On the sample distribution of the adjusted coefficient of determination \(TraditionalForm\`R\&_\^2\) in OLS

Thomas Colignatus, July 2 2006
http://www.dataweb.nl/~cool

Summary

In specification search, one invariably meets with the coefficient of determination (R^2). Minimizing the sum of squares of the errors is equivalent to maximizing (R^2). Including more explanatory variables increases the value of (R^2). Thus there is a temptation to include variables, even at random. The "adjusted (R^2)", denoted as \(TraditionalForm\`R\&_\^2\), is a first step for correction for inclusion of variables. Maximizing \(TraditionalForm\`R\&_\^2\)is equivalent to minimizing the estimated variance of the errors. The paper first clarifies the definition of the sample parameter (ρ^2) and then derives the densities of (R^2) and  \(TraditionalForm\`R\&_\^2\) and their expected values as depending upon (ρ^2). From this follows the confidence interval for (ρ^2), and it is also shown that \(TraditionalForm\`R\&_\^2\) is "almost" an unbiased estimator of (ρ^2). The formal deductions are supplemented with algebraic, numerical and graphical routines that allow a quick grasp of the relationships.

The following loads the relevant routines of The Economics Pack. (ResetAll is an application of CleanSlate.)

Needs["Economics`Pack`"]

ResetAll

Economics[Estimate]

SetOptions[Plot, TextStyle -> {FontSize -> 11}];

?AdjustedRSquared

The noncentral F-ratio distribution NoncentralFRatioDistribution[(Subscript[n, 1]), (Subscript[n, 2]), lambda] is the distribution of the ratio of a noncentral chi-square random variable with noncentrality parameter λ and (Subscript[n, 1]) degrees of freedom to a central chi-square random variable with (Subscript[n, 2]) degrees of freedom.

?NoncentralFRatioDistribution

dens[u1_] = 
((1 - u1)^(n2/2 - 1)*u1^(n1/2 - 1)*Gamma[(n1 + n2)/2]*
   Hypergeometric1F1Regularized[(n1 + n2)/2, n1/2, (u1*θ)/2])/
  (E^(θ/2)*Gamma[n2/2])

?RSquaredDensity

RSquaredDensity[Phi,  r2,  n,  
      k, φ]

?RSquaredExpectation

?AdjustedR2Density

?AdjustedR2Expectation

ptes1 = PlotExpectedR2s[15, 3, PlotRange -> {{0.4, 0.5}, {0.4, 0.6}}, 
    Draw -> {Thickness[0.015]}];

Abs solves little

eR2[n_, k_, rho2_] = RSquaredExpectation[Phi, Abs, n, k, Phi[rho2]]

eR2adj[n_, k_, rho2_] = AdjustedR2Expectation[Phi, Abs, n, k, Phi[rho2]]

ptes3  =  
        PlotExpectedR2s[30,  6,  First  ->  eR2, 
          Last  ->  eR2adj ];

Introduction

In specification search, one invariably meets with the coefficient of determination (R^2). Minimizing the sum of squares of the errors e'.e (see below for all notation) is equivalent to maximizing (R^2). Including more explanatory variables increases the value of (R^2). Thus there is a temptation to include variables, even at random. The "adjusted (R^2)", denoted as \(TraditionalForm\`R\&_\^2\), is a first step for correction for inclusion of variables. Maximizing \(TraditionalForm\`R\&_\^2\)is equivalent to minimizing the estimated variance of the error (s^2) = e'.e / (n - k) (Green (1993:245)). \(TraditionalForm\`R\&_\^2\) rises if variables are dropped that have a t-value less than 1 (Green (1993:193)). More sophisticated would be the use of the Amemiya and Akaike criteria (Green (1993:245) or Campos, Ericsson and Hendry (2005)). Currently we only consider \(TraditionalForm\`R\&_\^2\).

True significance always relates to the content of the analysis, but it can be useful to have information on statistical significance as well. If some search criterion has not been increased by a statistically significant amount, then the search doesn't seem successful. An expression of significance can thus be based upon the sample distribution of  (s^2) or  \(TraditionalForm\`R\&_\^2\). Below, the focus is on the sample distribution of \(TraditionalForm\`R\&_\^2\).

R.A. Fisher 1915 derived the distribution for the case when the explanatory variables are normally distributed (random causes). Currently, we are interested in the case of Ordinary Least Squares where the explanatory variables are given (fixed causes). Note that the latter case still has one random element, the error, but we are interested in the relation between the explanators and the dependent variable.

While \(TraditionalForm\`R\&_\^2\) is a sample statistic, the question arises whether there is a parallel for the population. Earlier, Barten (1962) suggested an expression (discussed by Theil (1971:179)), that later was revised by Press and Zellner (1978) to be (ρ^2)= B / (B + n (σ^2)), with  B = β'X'X.β . Defining θ = B / (σ^2), a useful expression is (ρ^2)= θ / (θ + n), or θ = n (ρ^2) / (1 - (ρ^2)). Another useful definition is φ = θ / n = (ρ^2) / (1 - (ρ^2)) or (ρ^2) = φ / (φ + 1). We take \(TraditionalForm\`ρ\&^\^2\) = \(TraditionalForm\`R\&_\^2\) as an estimator of (ρ^2) and can show that it is almost unbiased. Basically, we have an unbiased estimator (s^2) for (σ^2), an unbiased estimator b for β, giving \(TraditionalForm\`B\&^\) = b'X'X.b - k (s^2) as unbiased E[\(TraditionalForm\`B\&^\)] = B (see below), and thus it is easy to define \(TraditionalForm\`θ\&^\) =  \(TraditionalForm\`B\&^\) / (s^2) = (\(TraditionalForm\`B\&^\) / (σ^2)) / ( (s^2) / (σ^2)), which gives us an estimator \(TraditionalForm\`ρ\&^\^2\) which indeed can be shown to be just \(TraditionalForm\`R\&_\^2\). The question that remains to discuss is what the distribution of this \(TraditionalForm\`ρ\&^\^2\) is and how it relates to (ρ^2).

Since B is conditional to the sample X, the (sample, "parent") value (ρ^2) = (ρ^2)(X) also depends upon the sample, so that (ρ^2) can best be called a "sample parameter". There is only a true population value for some limit value of X that may be difficult to define. Though we consider X fixed, it might be a more relevant question how we selected the X, since this affects the value of (ρ^2)(X) and the distribution of  \(TraditionalForm\`ρ\&^\^2\) = \(TraditionalForm\`R\&_\^2\). That said, having some value \(TraditionalForm\`R\&_\^2\)(n, k) and its distribution, we still could determine the probability how much a new \(TraditionalForm\`R\&_\^2\)(m, κ ) differs from it (or vice versa).

Carapeto and Holt (2003) clarify that P[\(TraditionalForm\`R\&_\^2\) < \(TraditionalForm\`x\&_\)] = P[1 - n/(n-k) (1 - (R^2)) < 1 - n/(n-k) (1 - x)] = P[(R^2) < x], for any x in [0, 1] and its transformed \(TraditionalForm\`x\&_\), so that (R^2) and \(TraditionalForm\`R\&_\^2\) have the same distribution. Yet it appears that  \(TraditionalForm\`ρ\&^\^2\) = \(TraditionalForm\`R\&_\^2\) is a more useful estimator of the relevant sample parameter (ρ^2), and, to determine its properties, it is sensible to work directly with its distribution.

Finally, there is a curious notion in the literature that the definition of (ρ^2) would be arbitrary and that, hence, all ensueing deductions would have arbitrary value. Supposedly, one could design various other formats and transformations, such as e.g. (η^2) = θ / (θ + n - k), and there would no sensible or rational criterion about what would be the best representation of the "true parameter" behind the multiple correlation coefficient. This notion is curious, since the given definition of (ρ^2) depends upon the population parameters β and σ, and hence is a derived sample parameter that clearly is as objective as those population parameters that one already has accepted. Thus it is strange to accept β and σ on the one hand, and to reject them on the other hand. Next, it can be shown that E[y'.y] = B + n (σ^2) = (θ + n) (σ^2), so there is a rationale to take this as the base of a ratio; and then the given definition of (ρ^2) nicely captures the same intuition, already expressed in (R^2) and \(TraditionalForm\`R\&_\^2\), that we want to summarize the relative weight of β versus σ. Finally, the value of the definition of (ρ^2) is enhanced by the observation, shown below, that E[\(TraditionalForm\`R\&_\^2\) ] ≈ (ρ^2). Thus the definition of (ρ^2) seems as arbitrary as the definition for the true mean of a sample.

Derivations

The ordinary least squares (OLS) model is y = X β + ϵ = X b + e, where β are fixed population parameters, ϵ ~ N[0, (σ^2)I ] and b = (((X' . X))^(-1)).X'.y. The idempotent and symmetric projection is P = X.(((X' . X))^(-1)).X' so that \(TraditionalForm\`y\&^\  = \ P\ \ y\), while M = I - P so that e = M.y = M.ϵ  since M.X = 0. There are N observations, K explanatory variables. The variables can best be measured in deviation from their means, so that, with the model redefined, the total degrees of freedom are n = N -1, X has k = K - 1 columns, while the error has n - k degrees of freedom. The unbiased estimated variance then is (s^2)= e 'e / (n - k).

As we have y ~ N[Xβ,  (σ^2) I ] and e ~ N[0, (σ^2) M] we also have Xb ~ N[Xβ,  (σ^2) P], noting that X.(b - β)  = P ϵ.

We find b - β = (((X' . X))^(-1)).X' ϵ  so that Var[b] = E[ (b - β)'.(b - β) ] = (σ^2)(((X' . X))^(-1)) and b ~ N[β, (σ^2)(((X' . X))^(-1))]. If we take V = Diag[ (σ^2)(((X' . X))^(-1))] which takes only the diagonal elements, then (V^(-1/2)). (b - β ) ~ N[0, I ], then, following "Student", we can perform t-tests on the significant difference of the parameters from zero, with the F-test from R.A. Fisher providing a joint test. When β ≠ 0, the off-diagonal elements become important. The following statistic (R^2) provides joint information on those β ≠ 0.

Since X'e = 0 the total sum of squares SST = y'.y = b'X'X.b  + e'e = y'.P.y  + e'e =SSR + SSE, with SSR = b'X'X.b = y'.P.y as the sum of squares of the regression and SSE the sum of squares of the error (or the sum of squares of the residuals). All these sums of squares are nonnegative. The coefficient of determination is (R^2) = SSR / SST (i.e. the square of the multiple correlation coefficient).

The target value of SSR is B = β'X'X.β  and the "population" coefficient of determination is (ρ^2)= B / (B + n (σ^2)). Note (again) that B  and (ρ^2) are sample functions but restricted by the population parameters. It may be useful to distinguish any sample from a fixed observed value P, so that we can consider E[SSR | P] etcetera. Since the latter condition is obvious, it is not necessary to always mention it.

To determine the sample distribution of (R^2) we must investigate the distribution of its components, check their possible dependence and then consider their ratio.

On SST

We first observe that E[SST] = E[ y'y ] = E[(X β + ϵ)'.(X β + ϵ)] = (B + n (σ^2)) so that E[y'y  / (σ^2)] = θ + n = n (φ + 1) = n / (1 - (ρ^2)). Note the infinities when (σ^2) = 0 or (ρ^2) = 1.

We have y ~ N[Xβ, (σ^2)I ] and, using ξ = Xβ / σ, then y / σ ~ N[ξ, I] and we find y'y / (σ^2) ~ (χ^2)[n, ξ'ξ] which is the non-central chi-square with non-centrality parameter λ = ξ'ξ . The expectation of this distribution is n + λ and thus we verify that ξ'ξ = β'X'X.β / (σ^2) = θ. Hence y'y / (σ^2) ~ (χ^2)[n, θ], where θ = n (ρ^2) / (1 - (ρ^2)).

On SSR

Subsequently, SSR = b'X'X.b = (P ϵ  + Xβ)'. (P ϵ  + Xβ) = ϵ' P ϵ + 2 β'X' ϵ  + B.  The expected value of the SSR is E[SSR] = E[b'X'X.b] = E[ϵ' P ϵ ] + B, since E[ϵ] = 0. An idempotent matrix has only characteristics roots 1 or 0, hence the trace of P is equal to its rank, and thus E[b'X'X.b] = k (σ^2)  + B.

Using the unbiased (s^2), and using the independence of (s^2) and b (Green (1993:185)), we find unbiased estimator \(TraditionalForm\`B\&^\) = b'X'X.b - k (s^2), since E[\(TraditionalForm\`B\&^\)] = E[b'X'X.b - k (s^2)] = E[b'X'X.b] - E[k (s^2)] = (k (σ^2)  + B) -  k (σ^2)  =  B.

Thus, this other allocation of squares \(TraditionalForm\`B\&^\) =  b'X'X.b - k (s^2) = y'.y - n (s^2) gives us an unbiased estimator of B. This  \(TraditionalForm\`B\&^\)  can be called the "adjusted SSR" (SSRadj). Note that \(TraditionalForm\`B\&^\) / y'.y = (b'X'X.b - k (s^2)) / y'.y = (b'X'X.b / y'.y) - (k (s^2)  / y'.y) = (R^2) - k / (n - k) e'e / y'.y  = (R^2) - k / (n - k) (1 - (R^2)) = \(TraditionalForm\`R\&_\^2\), which expression is commonly called the adjusted (R^2).

This estimator \(TraditionalForm\`B\&^\) for B thus inspires the consideration of the estimator \(TraditionalForm\`ρ\&^\^2\) =  \(TraditionalForm\`B\&^\) / (\(TraditionalForm\`B\&^\) + n (s^2)) = (y'.y - n (s^2)) / y'.y = 1 - n (s^2) / y'.y = 1 - n/(n-k) SSE / SST = 1 - n/(n-k) (1 - (R^2)) = \(TraditionalForm\`R\&_\^2\), the adjusted (R^2) in other formula formats. Note that \(TraditionalForm\`ρ\&^\^2\) = \(TraditionalForm\`R\&_\^2\)= (1 - 1/Ftest) (R^2) where the Ftest-value = (n - k) / k . (R^2)/ (1 - (R^2)), which finds its use in testing the hypothesis that β = 0 or (R^2) = 0. Note also that  -k / (n - k) ≤ \(TraditionalForm\`R\&_\^2\) ≤ 1 and that \(TraditionalForm\`R\&_\^2\) > 0 iff (R^2) > k / n. Restated, (R^2) =  1 - (n - k) / n (1 - \(TraditionalForm\`R\&_\^2\)) and hence E[(R^2)] =  1 - (n - k) / n (1 -E[ \(TraditionalForm\`R\&_\^2\)]). Below we will show that E[\(TraditionalForm\`R\&_\^2\) ] ≈ (ρ^2) and hence E[(R^2)] ≈ 1 - (n - k) / n (1 - (ρ^2)). As we have 0 ≤ (ρ^2) ≤ 1 then it is logical that E[ (R^2)] ≥ k / n. We can already show that (lim_(θ  ->  0)) E[ (R^2)] = k / n. Note that SSR = b'X'X.b ≈ B + k (s^2) and certainly in expectation. Hence (R^2) ≈  (B + k (s^2)) /  (B + k (s^2)+ (n - k) (s^2)).  For B = 0 the latter collapses to k / n.  If β = 0 or (ρ^2) = 0 then the SSR  b'.X'.X.b of OLS can still generate a nonnegative value such that 0 < (R^2) < k / n, but we would take this as an indication that the true correlation is close to zero. Alternatively said, a condition for a successful regression (i.e. showing a contribution of the regressors) is that  b'.X'.X.b > k (s^2),  or \(TraditionalForm\`B\&^\)  >  0.

It follows that we actually are less interested in (R^2) and more interested in  \(TraditionalForm\`R\&_\^2\). Thus, our definition of (R^2) suggests an expression for a sample parameter (ρ^2), subsequently  \(TraditionalForm\`R\&_\^2\) appears to be an estimator of that sample parameter, and hence (R^2) drops from our considerations.

We would like E[ \(TraditionalForm\`ρ\&^\^2\) ] = (ρ^2). If it would be possible to do so, we would take the expectation of the numerator independently from the denominator, so that E[ \(TraditionalForm\`ρ\&^\^2\) ] = ? = E[\(TraditionalForm\`B\&^\) ] / E[\(TraditionalForm\`B\&^\) + n (s^2)] = B / (B + n (σ^2)) = (ρ^2). The terms are dependent however, and we cannot simply do that. But we can show below in a graphical manner that E[ \(TraditionalForm\`ρ\&^\^2\) ] ≈ (ρ^2), in fact very close.

As we have b'X'X.b - k (s^2) = y'.y - n (s^2) we can consider the distribution of either side of the equation.

Considering n (s^2), we know that  e'e / (σ^2) ~ (χ^2)[n - k] , so that n (s^2) is a n / (n - k) sum of those terms, so that n (s^2) / (σ^2) ~ n / (n - k) (χ^2)[n - k]  = (χ^2)[n] .

The ratio

To determine the sample distribution of \(TraditionalForm\`ρ\&^\^2\) = \(TraditionalForm\`R\&_\^2\) we consider the distribution of its components again.

Using  \(TraditionalForm\`ρ\&^\^2\) =  1 - n (s^2) / y'.y we find (with both n degrees of freedom):

\(TraditionalForm\`ρ\&^\^2\) ~   1 - (χ^2)[n]  / (χ^2)[n, θ], where  θ = n (ρ^2) / (1 - (ρ^2)).

However, these chisquares are dependent, and thus we cannot use known relations (such as on the F-distribution).

A successful strategy is to derive the distribution of \(TraditionalForm\`ρ\&^\^2\) = \(TraditionalForm\`R\&_\^2\) from that of (R^2), since  (R^2)= b'X'X.b / (b'X'X.b  + e'e) and we know that b and e are independent.

We have Xb ~ N[Xβ,  (σ^2) P] and let P C = C Λ, where (C^(-1)) = C '  are orthogonal eigenvectors, and where the diagonal matrix Λ has k non-zero eigenvalues.  In fact, since P is idempotent there are k values 1. Use z = C.X.b / σ and ζ = C.X.β / σ so that z ~ N[ζ, Λ ]. Thus z'z ~  (χ^2)[k, ζ '. ζ]. Obviously z'z = (b'.X' C'.).C.X.b / (σ^2) = b'X'X.b / (σ^2) and ζ '. ζ = B / (σ^2) = θ, so that b'X'X.b / (σ^2) ~  (χ^2)[k, θ]. We used the property of the eigenvalue decomposition, but we don't have to actually calculate them, since they disappear again.

Thus (R^2) = b'X'X.b / (b'X'X.b  + e'e) ~  (χ^2)[k, θ]  / ((χ^2)[k, θ]  + (χ^2)[n - k] ) where θ = n (ρ^2) / (1 - (ρ^2)).

Thus (R^2) = f / (f + 1) = 1 / (1 + 1 / f ) where f = (R^2) / (1 - (R^2)) so that r = f (n - k) / k  ~ F[k, n - k,  θ] which is the noncentral F-distribution, with in the numerator a noncentral chi-square random variable with noncentrality parameter λ = θ and (Subscript["n", 1])= k degrees of freedom and in the denominator a central chi-square random variable with (Subscript["n", 2]) = n - k degrees of freedom. We actually have already seen that transformation, see φ, though φ relates to the sample parameter of adjusted  (R^2) while f is a sample statistic of (R^2) itself. The Jacobian then is ∂r / ∂ (R^2) = 1 / (((1-(R^2)))^2) (n - k) / k .

Derivation using Mathematica

In the following, we will use the labels "density" and "expectation" for expressions in terms of φ, and the labels "dens" and "exp" for expressions in terms of (ρ^2).

Using the noncentral F distribution

The following functions are useful:

?Rho2

?Phi

?Theta

Note that we might replace \(TraditionalForm\`θ\&^\)  by Floor[\(TraditionalForm\`θ\&^\)] so that some functions simplify. For high values of rho2 the error need not be large.

Plot[Theta[rho2, 1], {rho2, 0, .95}];

Mathematica: "The noncentral F-ratio distribution NoncentralFRatioDistribution[(Subscript["n", 1]), (Subscript["n", 2]), lambda] is the distribution of the ratio of a noncentral chi-square random variable with noncentrality parameter λ and (Subscript["n", 1]) degrees of freedom to a central chi-square random variable with (Subscript["n", 2]) degrees of freedom."

Note that the f-value that is submitted to the function in Mathematica must itself be corrected for the degrees of freedom (which causes double entry).

For us, the substitutions result into:

PDF[NoncentralFRatioDistribution[k, n - k, n φ], 
r2 / (1 - r2) * (n - k) / k] / (1 - r2)^2 * (n - k) / k

When we use simplification, then Mathematica increases complexity by creating factors like ((R^2)- 1)^n, which are, indeed, in the complex plane - which we are not interested in. For this reason we use a2 = 1 - r2. (PM. Note that this is only a simplification procedure in Mathematica and not a change of distribution that would require a Jacobian.)

res = PDF[NoncentralFRatioDistribution[k, n - k, n φ], 
r2 / a2 * (n - k) / k] / a2^2 * (n - k) / k

A first round of simplification:

res2 = Simplify[res  //  PowerExpand,  
        Assumptions -> (n > 2 && k > 0 && 
              k < n && ((n | k) ∈ 
                  Integers) && (0 ≤ r2 ≤ 
                  1) &&  (φ  ∈  
                  Reals)  &&   φ ≥ 0 && 
              Re[φ] ≥ 0 && 
              Im[φ] == 0  &&  
              0  ≤  a2  ≤  
                1)]

A second round of simplification:

res3 = Simplify[res2  //  PowerExpand,  
        Assumptions -> (n > 2 && k > 0 && 
              k < n && ((n | k) ∈ 
                  Integers) && (0 ≤ r2 ≤ 
                  1) && φ ≥ 0 && 
              Re[φ] ≥ 0 && 
              Im[φ] == 0  &&  
              0  ≤  a2  ≤  
                1)]

And now we can safely substitute 1 - r2 and find the much desired density, in an acceptable form.

res3  /.  a2  ->  1  -  r2

The Estimate` package has already defined this funtion.

RSquaredDensity[Phi,  r2, n, k, φ]

It can be verified that it is a density.

Integrate[RSquaredDensity[Phi, r2, n, k, φ], {r2, 0, 1}, Assumptions -> 
    (n > 2 && k > 0 && k < n && ((n | k) ∈ Integers) && (φ ∈ Reals) && φ >= 0)]

The density in terms of parameter (ρ^2) is:

RSquaredDensity[r2, n, k, ρ2]

?RSquaredDensity

For example, when (ρ^2)  =  50  %  then  φ  =  
          1.

p1 = Plot[RSquaredDensity[r2, 15, 3, .5], {r2, 0, 1}, AxesLabel -> {"\!\(R\^2\)", density}];

p2 = Plot[RSquaredDensity[r2, 15, 3, .9], {r2, 0, 1}, AxesLabel -> {"\!\(R\^2\)", density}, PlotRange -> All];

Show[p1, p2, PlotRange -> All];

For larger values of n the calculation becomes slower. This happens when we double the sample size to 30.

Plot[RSquaredDensity[r2, 30, 3, .9], {r2, 0, 1}, PlotRange -> All, AxesLabel -> {\!\("R"\^2\), density}];

Expected value

Adding the right conditions gives the expected value - which again contains complex terms but cannot be simplified.

ev = Integrate[r2 RSquaredDensity[Phi, r2, n, k, φ], {r2, 0, 1}, Assumptions -> 
    (φ ∈ Reals && r2 ∈ Reals && n > 2 && k > 0 && k < n && ((n | k) ∈ Integers) && (0 <= r2 <= 1) && φ >= 0 && Re[φ] >= 0 && Im[φ] == 0)]

ev2 = FullSimplify[ev, 
Assumptions -> (n > 2 && k > 0 && k < n && (n | k) ∈ Integers && φ ∈ Reals && Re[φ] >= 0)]

We are however not interested in the complex plane. Hence we define a general function that allows us to select the real parts, the imaginary part and the absolute value. The best shortcut is to take the real part.

RSquaredExpectation[Phi,  func,  n,  k,  φ]

RSquaredExpectation[n,  k,  ρ2]

The expected value of (R^2) has a lower limit k / n but we cannot find this formally.

Limit[RSquaredExpectation[Phi, Identity, n, k, φ], φ -> 0, Assumptions -> 
    (n > 2 && k > 0 && k < n && ((n | k) ∈ Integers) && (φ ∈ Reals) && Re[φ] >= 0)]

The formal Limit does not work but in numerical cases it does. Note that (ρ^2) -> 0 iff φ -> 0.

RSquaredExpectation[Phi, Identity, 15, 3, φ]

Limit[%, φ -> 0]

We recognise this as k / n, and remember that \(TraditionalForm\`R\&_\^2\) > 0 iff (R^2) > k / n.  When E[\(TraditionalForm\`R\&_\^2\) ] ≈ (ρ^2) and  0 ≤ (ρ^2) ≤ 1 then it is logical that E[ (R^2)] ≥ k / n.

Let us switch to ρ2.

RSquaredExpectation[n,  k,  ρ2]

This is a plot. The diagonal gives the true value, and E[(R^2)] deviates.

pexpr2 = Plot[{rho2, RSquaredExpectation[15, 3, rho2]}, {rho2, 0, 1}, AxesLabel -> {"\!\(ρ\^2\)", "Exp[\!\(R\^2\)], n = 15, k = 3"}];

If the plot really gives a line and not merely a close approximation of a line, then E[(R^2)] = (ρ^2) + k/n (1 - (ρ^2)) and the bias is k/n (1 - (ρ^2)). One wonders whether one couldn't find a simpler expression for the density such that 
(∫_0^1) (R^2) f(
            (R^2),  (ρ^2)) d(R^2) = (ρ^2) + k/n (1 - (ρ^2)), of course with the proper confidence intervals.

Confidence interval

If we can find the CDF then we can find the confidence interval for ρ2 given an observed value of R2 as CDF[0, R2 ; ρ2Lower] = 1 - 0.025 to CDF[0, R2 ; ρ2Higher] = 0.025. However, Mathematica fails the integral for the CDF.

Integrate[ 
      RSquaredDensity[Phi,  r2, n, k, φ],  {r2,  0,  x}, 
      Assumptions -> (n > 2 && k > 0 && 
            k < n && ((n | k) ∈ 
                Integers) && (φ ∈ 
                Reals) && φ ≥ 0)]

For a numerical case, suppose that for n = 15 and k = 3 we observe a value R2 = .75. Then we can find a numerical solution:

RSquaredCI[ .75,  15,  3,   .05]

Note that we better use the adjusted R2 to find the confidence interval. But we cannot simply substitute it in the above function since that function uses the distribution of R2. Thus it is erroneous to use:

RSquaredCI[AdjustedRSquared[ .75,  15,  3],  15,  
      3,   .05]

The adjusted (R^2)

As \(TraditionalForm\`ρ\&^\^2\) = \(TraditionalForm\`R\&_\^2\) = 1 - n / (n-k) (1 - (R^2)) is a simple linear transformation, so that (R^2) =  1 - (n - k) / n (1 - \(TraditionalForm\`R\&_\^2\)) and Jacobian ∂ (R^2) / ∂  \(TraditionalForm\`R\&_\^2\) = (n - k) / n, we directly find the density of \(TraditionalForm\`R\&_\^2\).

(n-k)/n * RSquaredDensity[Phi, 1 - (n-k)/n (1 - r2adj), n, k, φ]

The package has already implemented this in:

AdjustedR2Density[Phi,  r2adj,  n,  k,  φ]

AdjustedR2Density[r2adj,  n,  k,  ρ2]

This provides a numerical example for \(TraditionalForm\`ρ\&^\^2\) when (ρ^2) = 50%.

p3 = Plot[AdjustedR2Density[r2, 15, 3, .5], {r2, -3/(15-3), 1}, PlotRange -> All, AxesLabel -> {"\!\(R\_adj\%2\)", density}];

This provides a numerical example for \(TraditionalForm\`ρ\&^\^2\) when (ρ^2) = 90%.

p4 = Plot[AdjustedR2Density[r2, 15, 3, .9], {r2, -3/(15-3), 1}, PlotRange -> All, AxesLabel -> {"\!\(R\_adj\%2\)", density}];

The mode of \(TraditionalForm\`R\&_\^2\) indeed is closer to the true value (ρ^2) than the mode of the distribution of (R^2).

Show[p1, p2, p3, p4, AxesLabel -> {"\!\(R\^2\) \nand \n\!\(R\_adj\%2\)", density}, PlotRange -> All];

Expected value of adjusted (R^2)

We can transform the expected value of (R^2), as E[\(TraditionalForm\`R\&_\^2\)] = 1 - n / (n - k) (1 - E[(R^2)] ].

AdjustedR2Expectation[Phi,  func,  n,  k,  φ]

AdjustedR2Expectation[n,  k,  ρ2]

And again we cannot find a simpler solution. But a plot is instructive. The plot for the adjusted Rsquare appears to be the diagonal, indicating unbiasedness.

pexpr2adj = Plot[{AdjustedR2Expectation[15, 3, rho2]}, {rho2, 0, 1}, AxesLabel -> {"\!\(ρ\^2\)", "Exp[\!\(R\_adj\%2\)], n = 15, k = 3"}, AxesOrigin -> {0, 0}, PlotRange -> All (*, PlotStyle -> {Thickness[.015] } *) ];

And combining the two plots shows that E[ \(TraditionalForm\`ρ\&^\^2\) = \(TraditionalForm\`R\&_\^2\)] ≈ (ρ^2) which confirms our intuition (or hope) that \(TraditionalForm\`R\&_\^2\) would be an unbiased estimator.

Show[pexpr2, pexpr2adj, AxesLabel -> {"\!\(ρ\^2\)", 
"Exp[\!\(R\^2\)], Exp[\!\(R\_adj\%2\)], \n n = 15, k = 3"}];

This routine plots for arbitrary n and k.

PlotExpectedR2s[10, 6, Draw -> {Thickness[.015]} ];

PlotExpectedR2s[35, 6];

The following plots indicate convergence problems, where where it doesn't help to switch between Re and Identity.

PlotExpectedR2s[30, 6, PlotRange -> {{0, 0.3}, {0, 0.3}}, PlotPoints -> 100];

RSquaredExpectation[30,  6,  rho2]

RSquaredExpectation[Phi,  Identity,  30,  6,  Phi[rho2]]

Plot[%,  {rho2,  0,  0.2},  
        AxesOrigin  ->  {0,  0}];

This strikes one as a bit strange.

PlotExpectedR2s[30, 6, PlotRange -> {{0, 0.1}, {0, 0.1}}, PlotPoints -> 100];

Anyway, this close-up suggests that adjusted R2 is not quite unbiased.

PlotExpectedR2s[30, 6, PlotRange -> {{0.6, 0.7}, {0.6, 0.71}}, PlotPoints -> 100, Draw -> {Thickness[.015]}];

Clearly, the issue is not fully settled due to (a) the imaginary values, (b) numerical precision anyway. To settle the issue, we would like more formal statements about the closeness.

AdjustedR2Expectation[n, 
        k, (ρ^2)] ≈  (ρ^2)

The appendices contain various efforts to get rid of the complex terms, to start with; all failed.

Confidence interval

Earlier we considered the confidence interval for (ρ^2) assuming that (R^2) = .75 had been observed for n = 15 and k = 3. Using the distribution for adjusted (R^2) we can now compare the original confidence interval with the one that we get when using the adjusted (R^2) and its density.

RSquaredCI[ .75,  15,  3,   .05]

AdjustedRSquared[ .75,  15,  3]

AdjustedRSquaredCI[%,  15,  3,   .05]

which is rather consistent.

Substituting the expected value into the density

When we have an observed \(TraditionalForm\`ρ\&^\^2\) then a next question is what the most likely distribution is from which it comes. A shortcut is to accept \(TraditionalForm\`R\&_\^2\) as the almost unbiased estimator and plot densities assuming this value for (ρ^2). This is formally doubtful and numerically risky but it allows quick visual comparison and gives an impression of the stituation.

?SetR2Plots

SetR2Plots[0.6875,  15,  3]

RSquaredDensity[Plot];

AdjustedR2Density[Plot];

PM. In the following combination the horizontal axis should actually mention both statistics.

Show[%,  %%];

PM. Now we assume normal errors and then derive the sample distribution of \(TraditionalForm\`R\&_\^2\). We might also work around, assume a Beta distribution for some (r^2), and then derive the distribution of the errors. Perhaps this gives formulas that are more manageable and still acceptable. Our objective is to make probability statements on the goodness of fit, and not to impose, at all cost, normality on the error.

Mathematica routine to combine it all

We can combine various results of the above discussion into a single routine RSquared.

Options[RSquared]

?RSquared

Application of the routine

Arbitrary example

This gives some arbitrary y and equally arbitrary \(TraditionalForm\`y\&^\). The conditions for an OLS estimator are not fulfilled so that \(TraditionalForm\`y\&^\) is rescaled.

RSquared[ {5. , 3, 2, 4, 2}, {4.5, 3, 3. , 3.9, 1.5}, 2] // 
      MatrixForm

Above routine sets the density function. Any value of (ρ^2) can be entered but the default is \(TraditionalForm\`R\&_\^2\).

PM. There is an Indeterminate value at (R^2) = 0 so that we have to set the plot range.

RSquaredDensity[Plot,  
        PlotRange  ->  {0,  3.5}];

AdjustedR2Density[Plot,  
        PlotRange  ->  {{-0.2,  
              1},  {0,   3.5}}];

The formula for this particular density for  (ρ^2) = R2Adj is:

RSquaredDensity[Set,  R2]

The same example but 5 * y

If we consider 5 y but the same \(TraditionalForm\`y\&^\) then the routine adjusts the scale.

RSquared[5 *  {5. , 3, 2, 4, 2}, {4.5, 3, 3. , 3.9, 1.5}, 2] // 
      MatrixForm

The key properties are not affected.

Other results

The results of the fit are stored in a Results object, giving (y - \(TraditionalForm\`y\&_\)) and (\(TraditionalForm\`y\&^\) - \(TraditionalForm\`\(y\&^\)\&_\)).

Results[RSquared,  List]

Results[RSquared,  Fit]

Results[RSquared,  Sum]

Negative \(TraditionalForm\`R\&_\^2\)

Negative adjusted RSquareds are a common phenomenon.

RSquared[{5. , 3, 2, 4, 2}, {4.5, 4, 3.5, 2, 1}, 2] // 
      MatrixForm

Note that we do not reject the Null Hypothesis that the RHS doesn't provide an explanation.

RSquaredDensity[Plot,  
        PlotRange  ->  {0,  3.5}];

When there is no real fit

Ideally, and always when the sum of squares is minimized,  \(TraditionalForm\`y\&_\) = \(TraditionalForm\`\(y\&^\)\&_\) and (\(TraditionalForm\`y\&^\) - \(TraditionalForm\`\(y\&^\)\&_\)).e = 0. If Fit -> True (default) then the estimate is automatically rescaled, under the assumption that the problem is merely numerical precision and rounding errors (otherwise raise K by 1). If the option Fit -> False is used then the test on those conditions will indicate whether these hold.

RSquared[{5. , 3, 2, 4, 2}, {4.5, 4, 3.5, 2, 1}, 2,  
        Fit  ->  False] // MatrixForm

If one wishes, once can force the use of the raw input, neglecting the zero conditions. The R2 and R2Adj are based upon 1 - e'e / y'y (but have little theoretical meaning if the conditions are not satisfied).

RSquared[{5. , 3, 2, 4, 2}, {4.5, 4, 3.5, 2, 1}, 2,  
        Fit  ->  False,  ZeroTest  ->  False] // 
      MatrixForm

Conclusion

Above relations are not included in standard references such as Theil (1971) or Green (1993). For specification search it seems advisable that they are more standardly included. Of course, the discussion above makes some assumptions here, and these may not be suitable, but presently they seem consistent.

The advantage of the discussion above is not merely that these relations have been stated but also that they are available as routines, with numbers with known accuracy and with formulas and graphs (that are the (temporary) human advantage over digital machines).

Appendices

Graphical decomposition

The density of R2 consists of a beta distribution (depending upon n and k) times another term (depending upon rho2). It doesn't seem that one can easily simplify the expression.

PDF[BetaDistribution[a,  b],  x]

betaR2[r2_, n_, k_]  = 
      PDF[BetaDistribution[k/2,  (n - k)/2],   r2]

densityR2[r2, n, k, φ]

prt[r2_, n_, k_, φ_]  = 
      Hypergeometric1F1[n/2, 
          k/2, (n*r2*φ)/
            2]/(E^((n*φ)/2))

Plot[prt[r2,  15, 3,  1],  {r2,  0,  1},  
        PlotRange  ->  All];

Plot[betaR2[r2,  15, 3],  {r2,  0,  1}];

densityR2[r2,  15,  3,  1]

Plot[
        betaR2[r2,  15, 3] * prt[r2,  15, 3,  1],  {r2,  0,  
          1}];

Plot[prt[r2,  15, 3,  phi[ .9]],  {r2,  0,  1},  
        PlotRange  ->  All];

Another way to find the same result

The awkward result on the expectation causes us to look for alternative ways to derive our results. It appears that we find the same results, which adds to our confidence, but, the results remain awkward.

We had  (R^2) = b'X'X.b / (b'X'X.b  + e'e) ~  (χ^2)[k, θ]  / ((χ^2)[k, θ]  + (χ^2)[n - k] ) where θ = n (ρ^2) / (1 - (ρ^2)).

Subsequently we follow the same calculation scheme in MATH2011 (see the references). There, the property is derived, using (χ^2)[p]  = γ[2, p / 2], that for two variates y1 ~ γ[2, n1 / 2] and y2 ~ γ[2, n2 / 2] then y1 / (y1 + y2) ~ Beta[n1, n2]. For us, one of the distributions is a non-central chisquare, but the derivation scheme is the same. Thus we select u1 = y1 / (y1 + y2) and u2 = y1 + y2, so that y1 = u1 u2 and y2 = u2 (1 - u1) and the Jacobian is u2 (which is, given the squares, always nonnegative). Later we substitute n1 = k and n2 = n - k. A special condition is n2 ≤ n1 or n - k ≤ k or k ≥ n / 2.

?ChiSquareDistribution

?NoncentralChiSquareDistribution

PM Mathematica notes: "The sum of the squares of n normally distributed random variables with variance (σ^2) = 1 and nonzero means follows a noncentral chi-square distribution NoncentralChiSquareDistribution[n, lambda]. The noncentrality parameter λ is the sum of the squares of the means of the random variables in the sum. Note that in various places in the literature, λ/2 or (Sqrt[λ]) is used as the noncentrality parameter."

g1[u1_, u2_] = PDF[NoncentralChiSquareDistribution[n1, θ], u1 u2] *
PDF[ChiSquareDistribution[n2], u2 (1 - u1)] * u2

Where the symbol \(TraditionalForm\`\(\_0\)\(F\&~\_1\)\)is the Hypergeometric0F1Regularized function as implemented in Mathematica.

Integrating out u2, we have the density of (R^2).

dens[u1_] = Integrate[g1[u1, u2], {u2, 0, ∞}, Assumptions -> 
    (n1 > 2 && n2 > 2 &&  ((n1 | n2) ∈ Integers) && 0 <= u1 <= 1 && (θ ∈ Reals) &&
   Re[θ] >= 0 && Im[θ] == 0)]

Where the symbol \(TraditionalForm\`\(\_1\)\(F\&~\_1\)\)is the Hypergeometric1F1Regularized function as implemented in Mathematica.

We can check formally that the integral is 1.

Integrate[dens[u1], {u1, 0, 1}, Assumptions -> 
    (n1 > 2 && n2 > 2 && ((n1 | n2) ∈ Integers) && 0 <= u1 <= 1) && (Re[θ] >= 0)]

We even find the expected value.

expval[u1_] = Integrate[u1 dens[u1], {u1, 0, 1}, Assumptions -> 
    (n1 > 2 && n2 > 2 && ((n1 | n2) ∈ Integers) && (0 <= u1 <= 1) && (θ ∈ Reals) && (Re[θ] >= 0) && Im[θ] == 0)]

Change to (R^2)

Taking n1 = k and n2 = n - k, and using proper variable r2 to indicate the (R^2). Also use φ = θ / n.

g1[r2_, n_, k_, φ_] = g1[r2, u2] /. {n1 -> k, n2 -> n-k, θ -> n φ} // FullSimplify

Integrating with these variables gives the same result as when we would directly substitute them in the result above, which indicates that the domains of n1 and n2 might be less relevant.

The main effect of this alternative route is that the Beta expression now is expressed in Gamma's.

densityR2[r2_, n_, k_, φ_] = Integrate[g1[r2, n, k, φ] , {u2, 0, ∞}, Assumptions -> 
    (n > 2 && k > 0 && k < n && ((n | k) ∈ Integers) && 0 <= r2 <= 1 && (φ ∈ Reals) && (Re[φ] >= 0) && (Im[φ] == 0))]

If we take k > 0 the integral is 1.

Integrate[densityR2[r2, n, k, φ], {r2, 0, 1}, Assumptions -> 
    (n > 2 && k > 0 && k < n && ((n | k) ∈ Integers) && φ >= 0)]

And the expected value is not different.

Note on the use of Mathematica

Mathematica has been extremely helpful in deriving these relations (1) by being a database on known mathematical relations, (2) for symbolic manipulation, (3) by supporting quick numerical examples, with worriless accuracy, (4) by its graphical interface, (5) by providing speed and convenience, (6) by being flexible so that routines could be taylored. Needless to say, pen and paper proved indispensable as well, but those didn't provide the other points above.

The focus was on deriving these relations. Much work is not shown here. For example, at one stage, I used Mathematica to derive the distribution of N[μ, (σ^2)] / (N[μ, (σ^2)] + (χ^2)[n]) but later I realized that I hadn't properly thought about interdependence, that is, I assumed independence for the derivation, but the assumptions didn't fit the problem. One can wonder whether the good features of Mathematica provide a sink, i.e. that one stops thinking properly and proceeds with a line of thought since it is so easy (that is, relatively) to continue. However, once one reaches a result, it is standard to review all steps, and errors are identified. Alternatively, if dead alleys are not investigated, one might as well give up the whole endeavour. Thus in this case I think that it is fortunate that Mathematica provided the environment to identify an intermediate result, so that it was possible to confidently reject it, and to refocus on the main line of reasoning.

I did notice that Mathematica is not strong on integrals of standard statistical functions. For example the expectation of k x when x ~ N[μ, (σ^2)] makes one eventually abort the operation.

Also, "explanations" provided by Mathematica can be non-informative. For example the "explanation" of (tanh^(-1)) leaves a person nonplussed, while it would be possible to give a small explanation.

?ArcTanh

The following can be noted as well.

pdf = PDF[NoncentralFRatioDistribution[n1, n2, λ], x n2 / n1]

Domain[NoncentralFRatioDistribution[n1, n2, λ]]

Mathematica is curiously unable to show what this default function integrates to. As it is a PDF, one would expect the value 1.

Integrate[
      pdf,  {x,  
        0,  ∞},  
Assumptions -> (n1 > 
              0 && n2 > 
              2  && ((n1 | n2) ∈ 
                Integers) && (λ  ∈  
                Reals)  &&  (λ ≥ 
                0))]

Decomposition of the expected value

The following takes the different multiplicative components of the expected value, so that we can inspect their contribution.

ev = RSquaredExpectation[Phi, Identity, n, k, φ]

decomp = {2^(n/2 - 1),  E^((n*φ)/2), ((-n)*φ)^(n/2), (k - n)*Gamma[n/2] + 2*Gamma[n/2 + 1, -((n*φ)/2)] - k*Gamma[n/2, -((n*φ)/2)] }

decomprho2 = decomp /. {n -> 15, k -> 3, φ -> rho2}

A conclusion of these graphs below is that large values compensate for small values, and that the imaginary values cannot be neglected (though the Re[...] of the total remains the most sensible selection). PM. We don't need to plot the constant.

Table[Plot @@ {Re[decomprho2[[i]]], {rho2, 0, 1}, PlotRange -> All}, {i, 2, 4}];

We only plot the two complex expressions.

Table[Plot @@ {Im[decomprho2[[i]]], {rho2, 0, 1}, PlotRange -> All}, {i, 3, 4}];

Some more derivations on SSRadj

The following should help to clarify our intuitions about the distribution of  \(TraditionalForm\`ρ\&^\^2\).

To better understand the situation, we also consider the LHS of SSRadj = (b'X'X.b - k (s^2)) = (y'.y - n (s^2)).

((b' X' X . b)/((σ^2))) = \(TraditionalForm\`ϵ'\/σ\). P. (ϵ/σ) + 2 \(TraditionalForm\`β'\/σ\)X' . (ϵ/σ)  +  B / (σ^2)   with θ = B / (σ^2).

The part \(TraditionalForm\`ϵ'\/σ\). P. (ϵ/σ)  is a sum of squares of standard normal variables P ϵ / σ, and thus has the (χ^2)[k] distribution.

The part \(TraditionalForm\`β'\/σ\)X' . (ϵ/σ) is a linear sum (\(TraditionalForm\`β'\/σ\)X' ) of the standard normal distibution (ϵ/σ)~ N[0, I  ]. Thus  \(TraditionalForm\`β'\/σ\)X' . (ϵ/σ)  ~ N[0, θ ]. But we have two of those, thus  2\(TraditionalForm\`β'\/σ\)X' . (ϵ/σ)  ~ N[0, 2 θ ]. The constant can be included as the mean of this distribution, thus  2 \(TraditionalForm\`β'\/σ\)X' . (ϵ/σ)  +  B /  (σ^2)  ~ N[θ, 2 θ ]. Note that (χ^2)[p] has mean p and variance 2 p, so that N[θ, 2 θ ] approximates a (χ^2)[θ] distribution.

Thus ((b' X' X . b)/((σ^2))) ~ (χ^2)[k]  + N[θ, 2 θ ], in a shorthand notation as the sum of two random variates.

If we sum the distributions of the two variates we must check whether they are dependent or not. A sufficient condition for independence of the linear form  L. (ϵ/σ)  and the quadratic form \(TraditionalForm\`ϵ'\/σ\). P. (ϵ/σ) is that L P = 0. Substituting L = \(TraditionalForm\`β'\/σ\)X' we find L P = \(TraditionalForm\`β'\/σ\)X'  X.(((X' . X))^(-1)).X'   = L ≠  0, which indicates dependence (which is consistent with the notion that a sum of squares would always generate nonnegative values).

Denote the "independent remainder" as D[θ]. Thus, as we had some variables a ~ (χ^2)[k] and b ~ N[θ, 2 θ ], mutually dependent, we now consider the variables a ~ (χ^2)[k]  and  f ~  D[θ], now independent, such that a + b ⇔ a + f.  (Above we had the total as (χ^2)[k,  θ] , thus  f ~  D[θ] = (χ^2)[k,  θ]  - (χ^2)[k]  can be used.)

The numerator of  \(TraditionalForm\`ρ\&^\^2\) also contains a term (- k (s^2)) and we know that  e'e / (σ^2) ~ (χ^2)[n - k]

In summary, we find for the sample distribution of  \(TraditionalForm\`ρ\&^\^2\) =  \(TraditionalForm\`B\&^\) / (\(TraditionalForm\`B\&^\) + n (s^2)) that:

(a) numerator ((b' X' X . b)/((σ^2)))  - k (s^2)/ (σ^2)  ~ ( (χ^2)[k]  + D[θ] ) - (k/(n - k))(χ^2)[n - k]

~ (χ^2)[k]   - (k/(n - k))(χ^2)[n - k]  + D[θ]

~ (χ^2)[
        k  -  (k/(n - k)) (n - k)]   + D[θ]

~ D[θ]

(b) denominator  (((b' X' X . b)/((σ^2)))  - k (s^2)/ (σ^2) ) + (n (s^2)/ (σ^2)) = ((b' X' X . b)/((σ^2)))  + (n - k) (s^2)/ (σ^2) = (b'X'X.b + e'e) /  (σ^2)

(((b' X' X . b)/((σ^2)))  - k (s^2)/ (σ^2) ) + (n (s^2)/ (σ^2))  ~ D[θ, k] +  (n/(n - k)) (χ^2)[n - k]

~ D[θ] +  (χ^2)[n]

Hence  \(TraditionalForm\`ρ\&^\^2\)  ~ D[θ] /  (D[θ] +  (χ^2)[n] ).

Since (ρ^2) = θ / (θ + n) or θ = n (ρ^2) / (1 - (ρ^2)), thus D[θ] = Δ[(ρ^2), n]. Note that θ and Δ, though using Greek letters, are still sample functions, but restricted by population parameters.

Hence  \(TraditionalForm\`ρ\&^\^2\)  ~ Δ[(ρ^2), n] /  (Δ[(ρ^2), n] +  (χ^2)[n] ) where it may be noted that there is no explicit term for k.

We would like E[ \(TraditionalForm\`ρ\&^\^2\) ] = (ρ^2). If it is possible to so, which is the big IF, we select the f  above such that we can also take the expectation of the numerator independently from the denominator, so that E[ \(TraditionalForm\`ρ\&^\^2\) ] = E[\(TraditionalForm\`B\&^\) ] / E[\(TraditionalForm\`B\&^\) + n (s^2)] =  (ρ^2) and in particular such that D[θ] = Δ[(ρ^2), n] = (χ^2)[n]  (ρ^2) / (1 - (ρ^2)) = (χ^2)[n  (ρ^2) /(1-(ρ^2))]  = (χ^2)[θ] .

In that case  \(TraditionalForm\`ρ\&^\^2\)  ~ (χ^2)[θ] /  ((χ^2)[θ] +  (χ^2)[n] ).

Subsequently, we use the property that (χ^2)[p]  = γ[2, p / 2] and that for two variates x ~ γ[2, n / 2] and y ~ γ[2, m / 2] then x / (x + y) ~ Beta[n, m].

Thus:  \(TraditionalForm\`ρ\&^\^2\)  ~ Beta[θ, n ] . That is, if above assumptions on f  can be made. (The following shows that this runs into problems.)

An estimator is \(TraditionalForm\`θ\&^\) = n  \(TraditionalForm\`ρ\&^\^2\) / (1 -  \(TraditionalForm\`ρ\&^\^2\)), for convenience rounded down to an integer ("Floor[\(TraditionalForm\`θ\&^\) ]"), while negative \(TraditionalForm\`R\&_\^2\) can be set at a small positive value. It is a leap of faith to replace the sample parameter θ by its unbiased estimator \(TraditionalForm\`θ\&^\), in particular to describe the sample distribution of the \(TraditionalForm\`ρ\&^\^2\) that defines it, but the result is consistent.

Thus:  \(TraditionalForm\`ρ\&^\^2\)  ~ Beta[ \(TraditionalForm\`θ\&^\), n ] = Beta[Max[epsilon, Floor[n  \(TraditionalForm\`ρ\&^\^2\) / (1 -  \(TraditionalForm\`ρ\&^\^2\))]], n ].

Rather though, such a substitution causes an inaccuracy, and we rather have (re-introducing k):

Thus:  \(TraditionalForm\`ρ\&^\^2\)  ~ Beta[ \(TraditionalForm\`θ\&^\),  n-k ] = Beta[Max[epsilon, Floor[n  \(TraditionalForm\`ρ\&^\^2\) / (1 -  \(TraditionalForm\`ρ\&^\^2\))]],  n-k ].

But this still would not allow for negative values.

Thus the assumption that \(TraditionalForm\`ρ\&^\^2\)  ~ (χ^2)[θ] /  ((χ^2)[θ] +  (χ^2)[n] ) runs into problems, since \(TraditionalForm\`ρ\&^\^2\) can have negative values.

PM. We cannot use  \(TraditionalForm\`ρ\&^\^2\)  ~ (χ^2)[θ ] /  ((χ^2)[θ ] +  (χ^2)[n] ) = (χ^2)[θ ] / (χ^2)[θ + n ] and then use the F-distribution, since the terms in the ratio are dependent.

Aspects of the R2 density

On the hypergeometric function

?Hypergeometric0F1Regularized

PM Mathematica notes: "A limiting form of the confluent hypergeometric function which often appears is Hypergeometric0F1[a, z]. This function is obtained as the limit (Subscript[, 0])(Subscript[F, 1])(;a;z)=(lim_(q->∞))(Subscript[, 1])(Subscript[F, 1])(q;a;z/q). 
The (Subscript[, 0])(Subscript[F, 1]) function has the series expansion (Subscript[, 0])(Subscript[F, 1])(;a;z)=(∑_(k=0)^∞)1/(Subscript[(a), k]) (z^k)/k! and satisfies the differential equation zy''+ay'-y=0."

?Hypergeometric1F1Regularized

PM Mathematica notes: "The confluent hypergeometric function can be obtained from the series expansion (Subscript[, 1])(Subscript[F, 1])(a;b;z)=1+a  z/b+a(a+1)/b(b+1) (z^2)/2!+⋯=(∑_(k=0)^∞)(Subscript[(a), k])/(Subscript[(b), k]) (z^k)/k! . Some special results are obtained when a and b are both integers. If a<0, and either b>0 or b<a, the series yields a polynomial with a finite number of terms."  And: "The (Subscript[, 1])(Subscript[F, 1]) function can be written in the integral representation (Subscript[, 1])(Subscript[F, 1])(a;b;z)=Γ(b)/[Γ(b-a)Γ(a)](∫_( 0)^1)(e^(z  t))(t^(a-1)) (((1-t))^(b-a-1)) dt. " Also:

\(Hypergeometric1F1Regularized[a, b, 
        z] == ∑\+\(k = 0\)\%∞\( Pochhammer[a, k] z\^k\)\/\(\(\ \)\(Gamma[b + k] \(k!\)\)\)\)

Plot[Hypergeometric1F1Regularized[15/2, 13/2, 2 u1], {u1, 0, 1}];

Some notes on the gamma

ev2  =  
      RSquaredExpectation[Phi,  Identity,  n,  
        k,  φ]

Mathematica makes life complex. We are only interested in real values, but Mathematica makes sure that we don't forget about the complex plane, and they present us with generalizations which we might not need: "Gamma[z] is defined by the integral Γ(z)=(∫_0^∞)(t^(z-1))(e^(-t))dt. For positive integer n, Γ(n)=(n-1)! . Γ(z) can be viewed as a generalization of the factorial function, valid for complex arguments z." Also, Gamma[a, z] is defined by the integral Γ(a,z)=(∫_z^∞)(t^( a-1))(e^(-t))dt. Note that the z in Γ[z] determines the power whereas the a does in Γ[a, z], which doesn't look as consistent notation. Note that the incomplete Γ above uses negative values, and the documentation of Mathematica states: "Gamma[a, z] has a branch cut discontinuity in the complex z plane running from -∞ to 0". But we are only interested in real values, as our assumptions clearly stated that Im[θ] == 0.

MyGamma[z_,  a_]  :=  Gamma[a,  z]

For us, φ ≥ 0 and hence only the negative axis for -zz in Γ[z, -zz]  is relevant.

MyGamma[z_,  a_]  =  
      Gamma [z]  +  
        Integrate[ 
          t^(z - 1) E^(-t),  {t,  a,  0}]

Mathematica doesn't recognize this formally:

Gamma[z,  0]  ==  Gamma[z]

But for random numbers it does:

(rand  =  Random[];  {rand,  Gamma[rand,  0],  
        Gamma[rand],   Gamma[rand,  0]  ==  
          Gamma[rand]})

Note that this doesn't hold for the negative numbers

(rand  =  -Random[];  
{rand,  
        Gamma[rand,  0],  
        Gamma[rand],   Gamma[rand,  0]  ==  
          Gamma[rand]})

And neither for complex numbers.

(rand  =  I Random[];  {rand,  Gamma[rand,  0],  
        Gamma[rand],   Gamma[rand,  0]  ==  
          Gamma[rand]})

Consider the real part, the imaginary part, or the absolute value, for the negative axis.

Plot[Re[Gamma[3/2,  z]],  {z,  -10,  20}];

Plot[Abs[Gamma[3/2,  z]],  {z,  -10,  20}];

If we replace the incomplete Gamma under the assumption that the complex values below 0 would not be relevant, then the interaction with the other terms drops out, and we are left with a total value of 0.

ev2 /.  Gamma[z_,  zz_]  :>  Gamma[z]

%   //  FullSimplify

If we take only the Gamma terms:

h[n_, k_, φ_] = ((k - n)*Gamma[n/2] + 2*Gamma[n/2 + 1, -((n*φ)/2)] - 
    k*Gamma[n/2, -((n*φ)/2)])

h[15,  3,  φ]  //  FullSimplify

Plot[
        Re[h[15,  3,  φ]],  {φ,  0, 5}];

Plot[Abs[h[15,  3,  φ]],  {φ,  0,  .1},  
        PlotRange  ->  {0,   .1}];

Solve[h[n, k, φ]  ==  0,  φ]

Can Mathematica find the zero value ?

h[30,  6,  0]

SolveForZeroPhi[n_,  k_]  := 
Module[{obj},  
        obj  =  
          0  ==  
            Re[(((k - n)*Gamma[n/2] + 
                      2*Gamma[n/2 + 1, -((n*φ)/2)] - 
                      k*Gamma[
                          n/2, -((n*φ)/
                                2)]))]^2; 
Print[
          obj]; 
FindRoot  @@  {obj,  {φ,  .20,    .00,  1}}]

SolveForZeroPhi[30,  6]

This is no use. Various starting values give various results.

Some notes on the power of (- n φ)

Clear[f,  g]

f[θ_, n_] = (-θ)^(-n/2)

{ f[6,  6],  f[6,  7]}

Abs[%]  //  N

Just to be asure that we take account of the order.

g[θ_, 
        n_] = (((-θ)^n)^(1/
                2))^(-1)

{ g[6,  6],  g[6,  7]}

Re gives 0 for uneven x

Plot[Re[f[3,  x]],  {x,  0,  10}];

Im gives zero for even solutions

Plot[Im[f[2,  x]],  {x,  0,  10}];

Plot[Abs[f[2,  x]] ,  {x,  0,  10}];

Let us see how we can get rid of this superfluous complexity.

Series ?

It may be good to keep in mind that sometimes a series expansion clarifies matters. But not in this case, it would seem. Anyway, we can consider the series for small values of φ.

ev2  =  
      RSquaredExpectation[Phi,  Identity,  n,  
        k,  φ]

Series[ev2,  {φ,  0,  3}]  //  FullSimplify

%  /.  {n  ->  15,  k  ->  3}

which neither shows a way to get rid of the complex values.

References

Texts with a * were not available to the author. Texts with ** were too dense. Texts with *** are kept in mind. Note that some papers concern normally distributed regressors while we are interested in fixed regressors.

Anderson, T. W.  (1996), "R. A. Fisher and Multivariate Analysis",  Statistical Science, Vol. 11, No. 1, 20-34, http://projecteuclid.org/Dienst/UI/1.0/Summarize/euclid.ss/1032209662

* Barten, A.P. (1962), "Note on the unbiased estimation of the squared multiple correlation coefficient", Statistisca Neerlandica, 16, 151-163

*** Campos, J. N. R. Ericsson, and D. F. Hendry (2005), "General-to-specific Modeling:
An Overview and Selected Bibliography", Board of Governors of the Federal Reserve System
International Finance Discussion Papers Number 838, specifically section 4.5 on model selection criteria, http://www.nuff.ox.ac.uk/users/hendry/ElgarCEH.pdf

Carapeto, M. and W. Holt (2003), "On the equivalence of the distributions of R2 and adjusted R2", http://www.cass.city.ac.uk/faculty/mcarapeto/papers/R2.pdf

* Cramer, J.S. (1987), "Mean and variance of R2 in small and moderate samples", J. Econometrics 35,
253-266.

** Fisher, R. A. (1928), "The general sampling distribution of the multiple correlation coefcient". Proc. Roy. Soc. London Ser. A 121 654-673. See R.A. Fisher digital archive at
http://digital.library.adelaide.edu.au/coll/special/fisher/

Green, W.H. (1993), "Econometric analysis", Prentice Hall

Lohninger, H. (1999), "Teach/Me data analysis", http://www.vias.org/tmdatanaleng/cc_test_r.html (stochastic regressors)

MATH2011 Statistical distribution theory, Notes 8: Beta, F- and t-distribution,
http://www.maths.soton.ac.uk/teaching/units/math2011/handouts/Notes8.pdf

* Ohtani, K. (1994), "The density functions of R2 and, and their risk performance under asymmetric loss in misspecified linear regression models", Economic Modelling, 1994, vol. 11, issue 4, pages 463-471 
http://www.sciencedirect.com/science/article/B6VB1-46KC3HK-3/2/79f6781f3559cf10c775fa59ed68cabb

* Ohtani, K. (2000), "Bootstrapping R2 and adjusted R2 in regression analysis", Economic Modelling, 2000, vol. 17, issue 4, pages 473-483 http://econpapers.repec.org/scripts/search.asp?ft=Ohtani
see also http://www.econ.kobe-u.ac.jp/katudou/publication/paper/paper78-00.html

Ohtani, K. and H. Tanizaki (2004), "Exact distribution of (R^2)and adjusted (R^2)in a linear regression model with multivariate t error terms", J. Japan Statist. Soc. Vol 34 No 1 p101-109,
http://www.scipress.org/journals/jjss/pdf/3401/34010101.pdf

* Press, S.J. and A. Zellner (1978), "Posterior distribution for the multiple correlation coefficient with fixed regressors", Journal of Econometrics, 8, 307-321

* Steiger, http://www.interchange.ubc.ca/steiger/aboutr2.htm

Theil, H. (1971), "Principles of econometrics", Wiley

* Weisstein, E., http://www.mathworld.wolfram.com