24
$\begingroup$

What is the distribution of the coefficient of determination, or R squared, $R^2$, in linear univariate multiple regression under the null hypothesis $H_0:\beta=0$?

How does it depend on the number of predictors $k$ and number of samples $n>k$? Is there a closed-form expression for the mode of this distribution?

In particular, I have a feeling that for simple regression (with one predictor $x$) this distribution has mode at zero, but for multiple regression the mode is at a non-zero positive value. If this is indeed true, is there an intuitive explanation of this "phase transition"?


Update

As @Alecos showed below, the distribution indeed peaks at zero when $k=2$ and $k=3$ and not at zero when $k>3$. I feel that there should be a geometrical view on this phase transition. Consider geometrical view of OLS: $\mathbf y$ is a vector in $\mathbb R^n$, $\mathbf X$ defines a $k$-dimensional subspace there. OLS amounts to projecting $\mathbf y$ onto this subspace, and $R^2$ is squared cosine of the angle between $\mathbf y$ and its projection $\hat{\mathbf y}$.

Now, from @Alecos's answer it follows that if all vectors are random, then the probability distribution of this angle will peak at $90^\circ$ for $k=2$ and $k=3$, but will have a mode at some other value $<90^\circ$ for $k>3$. Why?!


Update 2: I am accepting @Alecos'es answer, but still have a feeling that I am missing some important insight here. If anybody ever suggests any other (geometrical or not) view on this phenomenon that would make it "obvious", I will be happy to offer a bounty.

$\endgroup$
  • 1
    $\begingroup$ Are you willing to assume error normality? $\endgroup$ – Dimitriy V. Masterov Dec 23 '14 at 1:17
  • 1
    $\begingroup$ Yes, I guess one has to assume it to make this question answerable (?). $\endgroup$ – amoeba Dec 23 '14 at 1:20
  • 1
    $\begingroup$ Have you checked this davegiles.blogspot.jp/2013/05/good-old-r-squared.html ? $\endgroup$ – Khashaa Dec 23 '14 at 2:35
  • 1
    $\begingroup$ @Khashaa: in fact, I have to admit that I did find that blogspot page before posting my question here. Honestly, I still wanted to have a discussion of this phenomenon on our forum, so pretended I did not see that. $\endgroup$ – amoeba Dec 23 '14 at 9:46
  • 1
    $\begingroup$ Strongly related CV question stats.stackexchange.com/questions/123651/… $\endgroup$ – Alecos Papadopoulos Dec 23 '14 at 13:14
31
$\begingroup$

For the specific hypothesis (that all regressor coefficients are zero, not including the constant term, which is not examined in this test) and under normality, we know (see eg Maddala 2001, p. 155, but note that there, $k$ counts the regressors without the constant term, so the expression looks a bit different) that the statistic

$$F = \frac {n-k}{k-1}\frac {R^2}{1-R^2}$$ is distributed as a central $F(k-1, n-k)$ random variable.

Note that although we do not test the constant term, $k$ counts it also.

Moving things around,

$$(k-1)F - (k-1)FR^2 = (n-k)R^2 \Rightarrow (k-1)F = R^2\big[(n-k) + (k-1)F\big]$$

$$\Rightarrow R^2 = \frac {(k-1)F}{(n-k) + (k-1)F}$$

But the right hand side is distributed as a Beta distribution, specifically

$$R^2 \sim Beta\left (\frac {k-1}{2}, \frac {n-k}{2}\right)$$

The mode of this distribution is

$$\text{mode}R^2 = \frac {\frac {k-1}{2}-1}{\frac {k-1}{2}+ \frac {n-k}{2}-2} =\frac {k-3}{n-5} $$

FINITE & UNIQUE MODE
From the above relation we can infer that for the distribution to have a unique and finite mode we must have

$$k\geq 3, n >5 $$

This is consistent with the general requirement for a Beta distribution, which is

$$\{\alpha >1 , \beta \geq 1\},\;\; \text {OR}\;\; \{\alpha \geq1 , \beta > 1\}$$

as one can infer from this CV thread or read here.
Note that if $\{\alpha =1 , \beta = 1\}$, we obtain the Uniform distribution, so all the density points are modes (finite but not unique). Which creates the question: Why, if $k=3, n=5$, $R^2$ is distributed as a $U(0,1)$?

IMPLICATIONS
Assume that you have $k=5$ regressors (including the constant), and $n=99$ observations. Pretty nice regression, no overfitting. Then

$$R^2\Big|_{\beta=0} \sim Beta\left (2, 47\right), \text{mode}R^2 = \frac 1{47} \approx 0.021$$

and density plot

enter image description here

Intuition please: this is the distribution of $R^2$ under the hypothesis that no regressor actually belongs to the regression. So a) the distribution is independent of the regressors, b) as the sample size increases its distribution is concentrated towards zero as the increased information swamps small-sample variability that may produce some "fit" but also c) as the number of irrelevant regressors increases for given sample size, the distribution concentrates towards $1$, and we have the "spurious fit" phenomenon.

But also, note how "easy" it is to reject the null hypothesis: in the particular example, for $R^2=0.13$ cumulative probability has already reached $0.99$, so an obtained $R^2>0.13$ will reject the null of "insignificant regression" at significance level $1$%.

ADDENDUM
To respond to the new issue regarding the mode of the $R^2$ distribution, I can offer the following line of thought (not geometrical), which links it to the "spurious fit" phenomenon: when we run least-squares on a data set, we essentially solve a system of $n$ linear equations with $k$ unknowns (the only difference from high-school math is that back then we called "known coefficients" what in linear regression we call "variables/regressors", "unknown x" what we now call "unknown coefficients", and "constant terms" what we know call "dependent variable"). As long as $k<n$ the system is over-identified and there is no exact solution, only approximate -and the difference emerges as "unexplained variance of the dependent variable", which is captured by $1-R^2$. If $k=n$ the system has one exact solution (assuming linear independence). In between, as we increase the number of $k$, we reduce the "degree of overidentification" of the system and we "move towards" the single exact solution. Under this view, it makes sense why $R^2$ increases spuriously with the addition of irrelevant regressions, and consequently, why its mode moves gradually towards $1$, as $k$ increases for given $n$.

$\endgroup$
  • 1
    $\begingroup$ Its mathematical. For $k=2$ the first parameter of the beta distribution (the "$\alpha$" in standard notation) becomes smaller than unity. In that case the Beta distribution has no finite mode, play around with keisan.casio.com/exec/system/1180573226 to see how the shapes change. $\endgroup$ – Alecos Papadopoulos Dec 23 '14 at 9:53
  • 1
    $\begingroup$ @Alecos Excellent answer! (+1) Can I strongly suggest that you add to your answer the requirement for the mode to exist? This is usually stated as $\alpha>1$ and $\beta>1$ but more subtly, it's ok if equality holds in one of the two ... I think for our purposes this becomes $k \geq 3$ and $n \geq k + 2$ and at least one of these inequalities is strict. $\endgroup$ – Silverfish Dec 23 '14 at 13:25
  • 2
    $\begingroup$ @Khashaa Except if theory demands it, I never exclude the intercept from the regression -it is the average level of the dependent variable, regressors or no regressors (and this level is usually positive, so it would be a foolishly self-created misspecification to omit it). But I always exclude it from the F-test of the regression, since what I care about is not whether the dependent variable has a non-zero unconditional mean, but whether the regressors have any explanatory power as regards deviations from this mean. $\endgroup$ – Alecos Papadopoulos Dec 24 '14 at 3:59
  • 1
    $\begingroup$ +1! Are there results for the distribution of $R^2$ for nonzero $\beta_j$? $\endgroup$ – Christoph Hanck Jan 6 '16 at 11:44
  • 1
    $\begingroup$ @ChristophHanck See also davegiles.blogspot.jp/2013/05/good-old-r-squared.html $\endgroup$ – Alecos Papadopoulos Jan 6 '16 at 14:19
18
$\begingroup$

I won't rederive the $\mathrm{Beta}(\frac{k-1}{2}, \, \frac{n-k}{2})$ distribution in @Alecos's excellent answer (it's a standard result, see here for another nice discussion) but I want to fill in more details about the consequences! Firstly, what does the null distribution of $R^2$ look like for a range of values of $n$ and $k$? The graph in @Alecos's answer is quite representative of what occurs in practical multiple regressions, but sometimes insight is gleaned more easily from smaller cases. I've included the mean, mode (where it exists) and standard deviation. The graph/table deserves a good eyeball: best viewed at full-size. I could have included less facets but the pattern would have been less clear; I have appended R code so that readers can experiment with different subsets of $n$ and $k$.

Distribution of R2 for small sample sizes

Values of shape parameters

The graph's colour scheme indicates whether each shape parameter is less than one (red), equal to one (blue), or more than one (green). The left-hand side shows the value of $\alpha$ while $\beta$ is on the right. Since $\alpha = \frac{k-1}{2}$, its value increases in arithmetic progression by a common difference of $\frac{1}{2}$ as we move right from column to column (add a regressor to our model) whereas, for fixed $n$, $\beta = \frac{n-k}{2}$ decreases by $\frac{1}{2}$. The total $\alpha + \beta = \frac{n-1}{2}$ is fixed for each row (for a given sample size). If instead we fix $k$ and move down the column (increase sample size by 1), then $\alpha$ stays constant and $\beta$ increases by $\frac{1}{2}$. In regression terms, $\alpha$ is half the number of regressors included in the model, and $\beta$ is half the residual degrees of freedom. To determine the shape of the distribution we are particularly interested in where $\alpha$ or $\beta$ equal one.

The algebra is straightforward for $\alpha$: we have $\frac{k-1}{2}=1$ so $k=3$. This is indeed the only column of the facet plot that's filled blue on the left. Similarly $\alpha < 1$ for $k<3$ (the $k=2$ column is red on the left) and $\alpha > 1$ for $k>3$ (from the $k=4$ column onwards, the left side is green).

For $\beta=1$ we have $\frac{n-k}{2}=1$ hence $k=n-2$. Note how these cases (marked with a blue right-hand side) cut a diagonal line across the facet plot. For $\beta > 1$ we obtain $k < n - 2$ (the graphs with a green left side lie to the left of the diagonal line). For $\beta < 1$ we need $k > n - 2$, which involves only the right-most cases on my graph: at $n=k$ we have $\beta=0$ and the distribution is degenerate, but $n=k-1$ where $\beta = \frac{1}{2}$ is plotted (right side in red).

Since the PDF is $f(x;\,\alpha,\,\beta) \propto x^{\alpha-1} (1-x)^{\beta-1}$, it is clear that if (and only if) $\alpha<1$ then $f(x) \to \infty$ as $x \to 0$. We can see this in the graph: when the left side is shaded red, observe the behaviour at 0. Similarly when $\beta<1$ then $f(x) \to \infty$ as $x \to 1$. Look where the right side is red!

Symmetries

One of the most eye-catching features of the graph is the level of symmetry, but when the Beta distribution is involved, this shouldn't be surprising!

The Beta distribution itself is symmetric if $\alpha = \beta$. For us this occurs if $n = 2k-1$ which correctly identifies the panels $(k=2, n=3)$, $(k=3, n=5)$, $(k=4, n=7)$ and $(k=5, n=9)$. The extent to which the distribution is symmetric across $R^2 = 0.5$ depends on how many regressor variables we include in the model for that sample size. If $k = \frac{n+1}{2}$ the distribution of $R^2$ is perfectly symmetric about 0.5; if we include fewer variables than that it becomes increasingly asymmetric and the bulk of the probability mass shifts closer to $R^2 = 0$; if we include more variables then it shifts closer to $R^2 = 1$. Remember that $k$ includes the intercept in its count, and that we are working under the null, so the regressor variables should have coefficient zero in the correctly specified model.

There is also an obviously symmetry between distributions for any given $n$, i.e. any row in the facet grid. For example, compare $(k=3, n=9)$ with $(k=7, n=9)$. What's causing this? Recall that the distribution of $\mathrm{Beta}(\alpha, \beta)$ is the mirror image of $\mathrm{Beta}(\beta, \alpha)$ across $x=0.5$. Now we had $\alpha_{k,n} = \frac{k-1}{2}$ and $\beta_{k,n} = \frac{n-k}{2}$. Consider $k'=n-k+1$ and we find:

$$\alpha_{k',n} = \frac{(n-k+1)-1}{2} = \frac{n-k}{2} = \beta_{k,n}$$ $$\beta_{k',n} = \frac{n-(n-k+1)}{2} = \frac{k-1}{2} = \alpha_{k,n}$$

So this explains the symmetry as we vary the number of regressors in the model for a fixed sample size. It also explains the distributions that are themselves symmetric as a special case: for them, $k' = k$ so they are obliged to be symmetric with themselves!

This tells us something we might not have guessed about multiple regression: for a given sample size $n$, and assuming no regressors have a genuine relationship with $Y$, the $R^2$ for a model using $k-1$ regressors plus an intercept has the same distribution as $1 - R^2$ does for a model with $k-1$ residual degrees of freedom remaining.

Special distributions

When $k=n$ we have $\beta=0$, which isn't a valid parameter. However, as $\beta \to 0$ the distribution becomes degenerate with a spike such that $\mathsf{P}(R^2 = 1)=1$. This is consistent with what we know about a model with as many parameters as data points - it achieves perfect fit. I haven't drawn the degenerate distribution on my graph but did include the mean, mode and standard deviation.

When $k=2$ and $n=3$ we obtain $\mathrm{Beta}(\frac{1}{2}, \, \frac{1}{2})$ which is the arcsine distribution. This is symmetric (since $\alpha = \beta$) and bimodal (0 and 1). Since this is the only case where both $\alpha < 1$ and $\beta < 1$ (marked red on both sides), it is our only distribution which goes to infinity at both ends of the support.

The $\mathrm{Beta}(1, \, 1)$ distribution is the only Beta distribution to be rectangular (uniform). All values of $R^2$ from 0 to 1 are equally likely. The only combination of $k$ and $n$ for which $\alpha = \beta =1$ occurs is $k=3$ and $n=5$ (marked blue on both sides).

The previous special cases are of limited applicability but the case $\alpha > 1$ and $\beta=1$ (green on left, blue on right) is important. Now $f(x;\,\alpha,\,\beta) \propto x^{\alpha-1} (1-x)^{\beta-1} = x^{\alpha-1}$ so we have a power-law distribution on [0, 1]. Of course it's unlikely we'd perform a regression with $k=n-2$ and $k>3$, which is when this situation occurs. But by the previous symmetry argument, or some trivial algebra on the PDF, when $k=3$ and $n > 5$, which is the frequent procedure of multiple regression with two regressors and an intercept on a non-trivial sample size, $R^2$ will follow a reflected power law distribution on [0, 1] under $H_0$. This corresponds to $\alpha=1$ and $\beta>1$ so is marked blue on left, green on right.

You may also have noticed the triangular distributions at $(k=5,n=7)$ and its reflection $(k=3,n=7)$. We can recognise from their $\alpha$ and $\beta$ that these are just special cases of the power-law and reflected power-law distributions where the power is $2-1=1$.

Mode

If $\alpha>1$ and $\beta>1$, all green in the plot, $f(x; \, \alpha, \, \beta)$ is concave with $f(0)=f(1)=0$, and the Beta distribution has a unique mode $\frac{\alpha-1}{\alpha+\beta-2}$. Putting these in terms of $k$ and $n$, the condition becomes $k>3$ and $n>k+2$ while the mode is $\frac{k-3}{n-5}$.

All other cases have been dealt with above. If we relax the inequality to allow $\beta=1$, then we include the (green-blue) power-law distributions with $k=n-2$ and $k>3$ (equivalently, $n>5$). These cases clearly have mode 1, which actually agrees with the previous formula since $\frac{(n-2)-3}{n-5}=1$. If instead we allowed $\alpha=1$ but still demanded $\beta>1$, we'd find the (blue-green) reflected power-law distributions with $k=3$ and $n>5$. Their mode is 0, which agrees with $\frac{3-3}{n-5}=0$. However, if we relaxed both inequalities simultaneously to allow $\alpha=\beta=1$, we'd find the (all blue) uniform distribution with $k=3$ and $n=5$, which does not have a unique mode. Moreover the previous formula can't be applied in this case, since it would return the indeterminate form $\frac{3-3}{5-5}=\frac{0}{0}$.

When $n=k$ we get a degenerate distribution with mode 1. When $\beta < 1$ (in regression terms, $n=k-1$ so there is only one residual degree of freedom) then $f(x) \to \infty$ as $x \to 1$, and when $\alpha < 1$ (in regression terms, $k=2$ so a simple linear model with intercept and one regressor) then $f(x) \to \infty$ as $x \to 0$. These would be unique modes except in the unusual case where $k=2$ and $n=3$ (fitting a simple linear model to three points) which is bimodal at 0 and 1.

Mean

The question asked about the mode, but the mean of $R^2$ under the null is also interesting - it has the remarkably simple form $\frac{k-1}{n-1}$. For a fixed sample size it increases in arithmetic progression as more regressors are added to the model, until the mean value is 1 when $k=n$. The mean of a Beta distribution is $\frac{\alpha}{\alpha+\beta}$ so such an arithmetic progression was inevitable from our earlier observation that, for fixed $n$, the sum $\alpha+\beta$ is constant but $\alpha$ increases by 0.5 for each regressor added to the model.

$$\frac{\alpha}{\alpha+\beta} = \frac{(k-1)/2}{(k-1)/2 + (n-k)/2} = \frac{k-1}{n-1}$$

Code for plots

require(grid)
require(dplyr)

nlist <- 3:9 #change here which n to plot
klist <- 2:8 #change here which k to plot

totaln <- length(nlist)
totalk <- length(klist)

df <- data.frame(
    x = rep(seq(0, 1, length.out = 100), times = totaln * totalk),
    k = rep(klist, times = totaln, each = 100),
    n = rep(nlist, each = totalk * 100)
)

df <- mutate(df,
    kname = paste("k =", k),
    nname = paste("n =", n),
    a = (k-1)/2,
    b = (n-k)/2,
    density = dbeta(x, (k-1)/2, (n-k)/2),
    groupcol = ifelse(x < 0.5, 
        ifelse(a < 1, "below 1", ifelse(a ==1, "equals 1", "more than 1")),
        ifelse(b < 1, "below 1", ifelse(b ==1, "equals 1", "more than 1")))
)

g <- ggplot(df, aes(x, density)) +
    geom_line(size=0.8) + geom_area(aes(group=groupcol, fill=groupcol)) +
    scale_fill_brewer(palette="Set1") +
    facet_grid(nname ~ kname)  + 
    ylab("probability density") + theme_bw() + 
    labs(x = expression(R^{2}), fill = expression(alpha~(left)~beta~(right))) +
    theme(panel.margin = unit(0.6, "lines"), 
        legend.title=element_text(size=20),
        legend.text=element_text(size=20), 
        legend.background = element_rect(colour = "black"),
        legend.position = c(1, 1), legend.justification = c(1, 1))


df2 <- data.frame(
    k = rep(klist, times = totaln),
    n = rep(nlist, each = totalk),
    x = 0.5,
    ymean = 7.5,
    ymode = 5,
    ysd = 2.5
)

df2 <- mutate(df2,
    kname = paste("k =", k),
    nname = paste("n =", n),
    a = (k-1)/2,
    b = (n-k)/2,
    meanR2 = ifelse(k > n, NaN, a/(a+b)),
    modeR2 = ifelse((a>1 & b>=1) | (a>=1 & b>1), (a-1)/(a+b-2), 
        ifelse(a<1 & b>=1 & n>=k, 0, ifelse(a>=1 & b<1 & n>=k, 1, NaN))),
    sdR2 = ifelse(k > n, NaN, sqrt(a*b/((a+b)^2 * (a+b+1)))),
    meantext = ifelse(is.nan(meanR2), "", paste("Mean =", round(meanR2,3))),
    modetext = ifelse(is.nan(modeR2), "", paste("Mode =", round(modeR2,3))),
    sdtext = ifelse(is.nan(sdR2), "", paste("SD =", round(sdR2,3)))
)

g <- g + geom_text(data=df2, aes(x, ymean, label=meantext)) +
    geom_text(data=df2, aes(x, ymode, label=modetext)) +
    geom_text(data=df2, aes(x, ysd, label=sdtext))
print(g)
$\endgroup$
  • 1
    $\begingroup$ Really illuminating visualization. +1 $\endgroup$ – Khashaa Dec 24 '14 at 8:16
  • $\begingroup$ Great addition, +1, thanks. I noticed that you call $0$ a mode when the distribution goes to $+\infty$ when $x\to 0$ (and nowhere else) -- something @Alecos above (in the comments) did not want to do. I agree with you: it is convenient. $\endgroup$ – amoeba Dec 24 '14 at 10:29
  • 1
    $\begingroup$ @amoeba from the graphs we'd like to say "values around 0 are most likely" (or 1). But the answer of Alecos is also both self-consistent and consistent with many authorities (people differ on what to do about the 0 and 1 full stop, let alone whether they can count as a mode!). My approach to the mode differs from Alecos mostly because I use conditions on alpha and beta to determine where the formula is applicable, rather than taking my starting point as the formula and seeing which k and n give sensible answers. $\endgroup$ – Silverfish Dec 24 '14 at 10:37
  • 1
    $\begingroup$ (+1), this is a very meaty answer. By keeping $k$ too close to $n$ and both small, the question studies in detail, and so decisively, the case of really small samples with relatively too many and irrelevant regressors. $\endgroup$ – Alecos Papadopoulos Dec 24 '14 at 12:39
  • $\begingroup$ @amoeba You probably noticed that this answer furnishes an algebraic answer for why, for sufficiently large $n$, the mode of the distribution is 0 for $k=3$ but positive for $k>3$. Since $f(x) \propto x^{(k-3)/2}(1-x)^{(n-k-2)/2}$ then for $k=3$ we have $f(x) \propto (1-x)^{(n-5)/2}$ which will clearly have mode at 0 for $n>5$, whereas for $k=4$ we have $f(x) \propto x^{1/2}(1-x)^{(n-6)/2}$ whose maximum can be found by calculus to be the quoted mode formula. As $k$ increases, the power of $x$ rises by 0.5 each time. It's this $x^{\alpha-1}$ factor which makes $f(0)=0$ so kills the mode at 0 $\endgroup$ – Silverfish Dec 24 '14 at 14:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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