8
$\begingroup$

I am looking at some simple regression models using both R and the statsmodels package of Python. I've found that, when computing the coefficient of determination, statmodels uses the following formula for $R^2$: $$ R^2 = 1 - \frac{SSR}{TSS}~~~~~~(\text{centered}) $$ where $SSR$ is the sum of squared residuals, and $TSS$ is the total sum of squares of the model. ("Centered" means that the mean has been removed from the series.) However, the same calculation in R yields a different result for $R^2$. The reason is that R seems to be calculating $R^2$ as: $$ R^2 = 1 - \frac{SSR}{TSS}~~~(\text{uncentered}) $$ So, what gives? Presumably there's some reason to prefer one over the other in certain situations. I haven't been able to find any information online about the cases where one of the above formulae should be preferred.

Can someone please explain why one is better than the other?

$\endgroup$
5
  • 1
    $\begingroup$ What is the model you use in R ? Does it include an intercept ? The second formula is the correct one when there's no interecept. $\endgroup$ Sep 10, 2012 at 21:50
  • $\begingroup$ Just to clarify, the statsmodels documentation states that rsquared is for a model that includes a constant. statsmodels.sourceforge.net/devel/generated/… $\endgroup$
    – Josef
    Sep 19, 2012 at 0:04
  • $\begingroup$ I understood that the definition of $R^2$ is different, what I didn't understand is why. $\endgroup$
    – BenDundee
    Sep 19, 2012 at 0:45
  • $\begingroup$ Could you give us the code that you're using? Especially in R, there are many different ways to calculate a regression model. $\endgroup$
    – analystic
    Sep 23, 2012 at 4:36
  • $\begingroup$ R code: lm(y ~ x+0), then summary, then look at the R^2 value. $\endgroup$
    – BenDundee
    Oct 1, 2012 at 13:08

5 Answers 5

11
$\begingroup$

As Stephane hinted in the comment, it is the difference between model with or without intercept that matters.

For what it is worth, here is code from one of my packages:

## cf src/library/stats/R/lm.R and case with no weights and an intercept
f <- object$fitted.values
r <- object$residuals
mss <- if (object$intercept) sum((f - mean(f))^2) else sum(f^2)
rss <- sum(r^2)
r.squared <- mss/(mss + rss)

Residuals are centered by design, leaves the fitted values with need to be centered in the intercept and not otherwise.

$\endgroup$
2
  • 3
    $\begingroup$ +1, nice point. Along those lines, it's worth taking a look at @cardinal's excellent answer here: removal-of-statistically-significant-intercept-term-boosts-r2. $\endgroup$ Sep 11, 2012 at 1:27
  • 1
    $\begingroup$ The colors suck, and I can't tell if I checked your response. I've clicked it a random number of times (by now), if you're not getting appropriate credit at this point, comment here and I'll click it once more :) $\endgroup$
    – BenDundee
    Sep 12, 2012 at 16:04
11
$\begingroup$

Ok, I thought I'd follow up on this. I've been struggling with the answers here a bit, and have come to some better understanding of the problem. For posterity, I also think that a full explanation of _why there are two different forms of this equation for R^2 would be beneficial to anyone who stumbles upon this thread. I don't know if this is common knowledge, or what--no one seems to explain (possibly a lot of people just don't know, or possibly it's so basic that it's expected that people 'just know') WHY there are two forms for R^2. This includes several sets of lecture notes by professors at major universities: perhaps I'm just not looking in the right places.

The reason for the two different equations above comes from the fact that you're comparing the model against the null hypothesis. The null hypothesis is "there exists zero relationship between the dependent and independent variables". This means you're taking the slope to be zero. Another way to say this is that you're comparing the regression model you build to a nested model with one fewer parameters.

Now, suppose we have a set of data with one independent variable (x) and one dependent variable (y). We have two choices:

  • We choose to model the relationship between x and y with a two parameter linear model (i.e., $\hat{y}_i = a_0 + a_1 \hat{x}_i + \epsilon_i$). The null hypothesis is $\hat{y}_i = a_0 + \epsilon_i$, and $\bar{y} \neq 0$ in general. Thus the appropriate form of $R^2$ to use is: $$ R^2 = 1- \frac{\sum_i (y_i - \hat{y}_i)^2}{\sum_i (y_i - \bar{y})^2} $$
  • We choose to model the relationship between x and y with a one parameter linear model, namely $\hat{y}_i = a_1 \hat{x}_i + \epsilon_i$. The null hypothesis is that there is no relationship between x and y, thus the correct null hypothesis is $\hat{y}_i = \epsilon_i$. In other words, the null hypothesis is just white noise. Clearly, $\mathbb{E}(y) = 0$, thus the correct form of $R^2$ is $$ R^2 = 1- \frac{\sum_i (y_i - \hat{y}_i)^2}{\sum_i y_i^2} $$

A good way to think about this is the following: suppose the null hypothesis were (100%) correct, and there truly were no relationship between x and y. What would we expect? If anything's fair, the answer is "We expect $R^2=0$."

  • In the case where we choose a two-parameter model, we expect that $\bar{y} = \hat{y}_i = a_0$. If this isn't obvious, try drawing the picture with the model value under the null hypothesis $\hat{y}_i$, the data point as $y_i$, and the average $\bar{y}$. If the model is correct (i.e., number of data points -> infinity), then you should be able to see graphically that $\bar{y} = \hat{y}_i = a_0$, in the case where the null hypothesis is true.
  • Conversely, using the same picture as above, $\hat{y}_i = \bar{y} = 0$. There's a slight rub here, because you have to worry about how these things go to zero. L'Hopital will tell you that, in this case at least, $\lim 0/0 = 0$, and everything is ok.

You can see why you get funny things happening with $R^2$ (like negative values) if you use the wrong form of the equation. I noticed it first because the statsmodels package in Python does one thing, and R does something else: it pains me to say, but R is right and statsmodels is wrong. (Well, not really "pains"...)

I would love for some feedback on this intuition. I have only found one reference where this is explained explicitly. Please see this pdf file (download here), Section 5.3.6. Additionally, the other linked answer on stackexchange alludes to this fact, but the reasoning wasn't completely clear to me (no offense to the person who answered the question, it is a very well-written response, and I can be dense at times! ).

Again, please correct my reasoning in the comments, and I will amend the post until it is acceptable.

$\endgroup$
2
  • 2
    $\begingroup$ That's a nice essay :) $\endgroup$ Sep 18, 2012 at 1:16
  • 2
    $\begingroup$ I figured it cost me enough hair...there's no need for someone else to go bald over this problem :) $\endgroup$
    – BenDundee
    Sep 18, 2012 at 14:48
0
$\begingroup$

I guess that, as Stephane hinted in the comment, all the difference is in how you describe your model. If you describe the same model, the r squared will be the same in both cases.

I will post some python code to show that afterward, but first a word of caution: statsmodels, with the OLS function do not add automatically the intercept, while the R formula will, so this may be the origin of your difference. If this make you feel uncomfortable, try to use the new formula syntax.

With the following code oyu can verify that they bot gives you the same value (statsmodels is actually tested to give the same results as R)

# create the data for statsmodels
import statsmodels.formula.api as smf
import pandas as pd
from pylab import arange, randn
x = arange(20)
y = x * 0.3 + randn(20)
data = pd.DataFrame({'x':x, 'y':y})

# create the data for Rpy
import rpy2.robjects as robjects
import pandas.rpy.common as com
r = robjects.r
from rpy2.robjects import FloatVector
from rpy2.robjects.packages import importr
stats = importr('stats')
base = importr('base')
robjects.globalenv["x"] = FloatVector(x)
robjects.globalenv["y"] = FloatVector(y)

# model without intercept
lm0 = stats.lm("y ~ x + 0")
s = base.summary(lm0)
print s.rx2("r.squared")[0], smf.ols('y ~ x + 0', data).fit().rsquared
print s.rx2("adj.r.squared")[0], smf.ols('y ~ x + 0', data).fit().rsquared_adj

# model with intercept
lm0 = stats.lm("y ~ x + 1")
s = base.summary(lm0)
print s.rx2("r.squared")[0], smf.ols('y ~ x + 1', data).fit().rsquared
print s.rx2("adj.r.squared")[0], smf.ols('y ~ x + 1', data).fit().rsquared_adj
$\endgroup$
2
  • $\begingroup$ It was indeed the source of the difference. I talked a bit with one of the statsmodels devs, and I think the updated version (I can't remember version numbers, by now!) behaves as R does---i.e., it calculates the correct R^2, given the model. $\endgroup$
    – BenDundee
    Mar 14, 2013 at 20:45
  • $\begingroup$ Yes, I guess that at the time of the question they could actually behave differently. I just chosen to answer just to give an update on the problem. Statsmodels is IMHO a good python alternative for statistics and deserve to be promoted a little :) btw, compliments for your answer, it was really educational $\endgroup$ Mar 14, 2013 at 21:13
0
$\begingroup$

My understanding is:

  • The uncentered R2 explains the explanatory power of all regressors (including the constant regressor).
  • The centered R2 explains the explanatory power of all nonconstant regressors.
$\endgroup$
0
$\begingroup$

When there is an intercept in the regression model, one can prove that

$ \sum (y_i-\bar{y})^2 = \sum (\hat{y}_i - \bar{y})^2 + \sum(y_i - \hat{y}_i)^2 $

The proof utilizes the two OLS first order conditions. And this equality ensures that $R^2$ is always in [0,1].

When there is no intercept, the equation above do not hold. Instead, we have $\sum y_i^2 = \sum \hat{y}_i^2 + \sum (y_i - \hat{y}_i)^2$ The proof utilizes the only OLS first order condition for this homogenous model. And the uncentered $R^2$ is always in [0,1]

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.