On the coefficient of determination
I have often been asked by my students why the $R^2$ measure is not a measure of comparison between two models.
In general, as will be shown below, the $R^2$ value is simply a measure of how well the model that you’ve used fits the data that you have. In particular, when comparing nested models estimated by OLS on the same data, $R^2$ can never decrease when additional regressors are added.
Let’s start with a basic linear regression model:
\[Y = \mathbf{X}\beta + \varepsilon,\]
where $Y$ is an $n\times1$ vector, $\mathbf{X}$ is an $n\times k$ matrix of regressors, $\beta$ is a $k\times1$ vector of coefficients, and $\varepsilon$ is an $n\times1$ vector of disturbances.
We assume the usual conditions required for the OLS estimator to have its standard properties. In particular, if we are interested in interpreting the coefficients causally, we require appropriate exogeneity assumptions. Causality itself, however, is not one of the assumptions of the classical linear regression model.
When the relevant exogeneity and rank conditions hold, the least squares estimator for $\beta$, denoted $b$, is consistent. If $\mathbf{X}^{T}\mathbf{X}$ has an inverse, then $b$ is estimated via the least squares normal equations as follows:
\[b = (\mathbf{X}^{T}\mathbf{X})^{-1}\mathbf{X}^{T}Y.\]
Denote the vector of OLS residuals by
\[e = Y - \mathbf{X}b.\]
Substituting the expression for $b$, we have
\[\begin{aligned}
e
&=Y-\mathbf{X}(\mathbf{X}^{T}\mathbf{X})^{-1}\mathbf{X}^{T}Y\\
&=\left(\mathbf{I}
-\mathbf{X}(\mathbf{X}^{T}\mathbf{X})^{-1}\mathbf{X}^{T}\right)Y.
\end{aligned}\]
Define the residual-maker matrix
\[\mathbf{M}_{X}
=
\mathbf{I}
-\mathbf{X}(\mathbf{X}^{T}\mathbf{X})^{-1}\mathbf{X}^{T}.\]
Then
\[e=\mathbf{M}_{X}Y.\]
Now, suppose we model the above with the addition of another variable:
\[Y = \mathbf{X}d + Zc + a,\]
where $Z$ is an $n\times1$ vector, $c$ is a scalar, and $a$ is the vector of residuals from the expanded model.
For a given value of $c$, the OLS estimator of the coefficient on $\mathbf{X}$ is
\[d
=
(\mathbf{X}^{T}\mathbf{X})^{-1}
\mathbf{X}^{T}(Y-Zc).\]
Using the expression for $b$, this can be written as
\[d
=
b-(\mathbf{X}^{T}\mathbf{X})^{-1}\mathbf{X}^{T}Zc.\]
Hence,
\[\begin{aligned}
a
&=Y-\mathbf{X}d-Zc\\
&=\mathbf{M}_{X}(Y-Zc)\\
&=e-\mathbf{M}_{X}Zc.
\end{aligned}\]
The Frisch–Waugh–Lovell theorem tells us that the coefficient $c$ in the expanded regression can be obtained by first removing the linear projection of $Z$ on $\mathbf{X}$, and the linear projection of $Y$ on $\mathbf{X}$, and then regressing the latter residuals on the former.
More importantly for our present purpose, the expanded model contains the original model as a special case. If we set $c=0$, then
\[Y=\mathbf{X}d+a,\]
and the original OLS solution is recovered by setting $d=b$.
The expanded model therefore minimizes the residual sum of squares over a parameter space that contains the parameter space of the original model. Consequently, the residual sum of squares from the expanded model cannot be larger than that from the original model:
\[a^{T}a \leq e^{T}e.\]
Finally, we arrive at $R^{2}$, also known as the coefficient of determination.
Suppose that the regression contains an intercept. Define $M^{0}$ as the $n\times n$ centering matrix
\[M^{0}
=
\mathbf{I}
-
\frac{1}{n}\mathbf{1}\mathbf{1}^{T},\]
where $\mathbf{1}$ is an $n\times1$ vector of ones.
The matrix $M^{0}$ transforms a vector into deviations from its sample mean. In particular,
\[Y^{T}M^{0}Y
=
\sum_{i=1}^{n}(Y_i-\bar Y)^2,\]
which is the total sum of squares.
For the regression
\[Y=\mathbf{X}b+e,\]
where $\mathbf{X}$ contains an intercept, the usual decomposition of the total sum of squares gives
\[Y^{T}M^{0}Y
=
b^{T}\mathbf{X}^{T}M^{0}\mathbf{X}b
+
e^{T}e.\]
The coefficient of determination is therefore
\[R^{2}
=
\frac{b^{T}\mathbf{X}^{T}M^{0}\mathbf{X}b}
{Y^{T}M^{0}Y}
=
1-
\frac{e^{T}e}
{Y^{T}M^{0}Y}.\]
Thus, $R^{2}$ measures the proportion of the sample variation in the dependent variable, relative to its sample mean, that is accounted for by the fitted regression.
Now consider the expanded model. Its residual sum of squares is $a^{T}a$, so its coefficient of determination is
\[R'^{2}
=
1-
\frac{a^{T}a}
{Y^{T}M^{0}Y}.\]
Since
\[a^{T}a\leq e^{T}e,\]
it follows immediately that
\[R'^{2}\geq R^{2}.\]
That is, one can never decrease the ordinary $R^{2}$ by adding regressors to a nested OLS model, provided the same dependent variable and observations are used.
This is precisely why $R^{2}$ is not, by itself, a useful criterion for deciding whether one model is better than another. If we simply add more regressors to a model, $R^{2}$ will necessarily increase or remain unchanged, even when the additional regressors provide little substantive value.
This does not mean that adding regressors is necessarily harmful. It means only that ordinary $R^{2}$ does not impose any penalty for model complexity.
For this reason, when comparing models, one might instead consider adjusted $R^{2}$, information criteria such as AIC or BIC, statistical tests of restrictions, or out-of-sample predictive performance, depending on the purpose of the model.
In particular, if the objective is prediction rather than in-sample fit, an increase in $R^{2}$ need not imply an improvement in predictive performance.