Comparing Regression Models and Multicollinearity, statistics lesson

A regression with one predictor usually leaves some variation in the response unexplained. Adding a second predictor and refitting by least squares on the same data lowers the residual sum of squares, or at worst leaves it unchanged. Does that make the larger model better? This lesson compares nested models with the partial F test and with adjusted R-squared, explains why closely related predictors make individual coefficients unstable or impossible to separate, and shows what searching through many candidate models does to the conclusions drawn afterward.

Key terms

Nested models

Definition A model is nested in a larger model when every mean function the smaller model allows is also allowed by the larger one. For linear models of the same response this holds when the column space of the smaller design matrix lies inside the column space of the larger one, for example when the smaller model is the larger one with some coefficients set to zero. Write $\text{SSE}$ for the residual sum of squares of a least-squares fit. For nested least-squares fits of the same response on the same observations,

\text{SSE}_{\text{larger}} \le \text{SSE}_{\text{smaller}} ,

so with the same total sum of squares $\text{SST}$ in both, the larger model’s $R^2 = 1 - \text{SSE}/\text{SST}$ is at least as large as the smaller model’s.

In plain words The larger model can reproduce the smaller fit by setting its extra coefficients to zero, so on the data used for fitting its squared error is never larger. That improvement is guaranteed by algebra, which is why a drop in $\text{SSE}$ or a rise in $R^2$ is not, on its own, evidence that the added term describes anything real.

Example The straight-line model $\beta_0 + \beta_1 x$ is nested in the quadratic model $\beta_0 + \beta_1 x + \beta_2 x^2$ , which gives the line when $\beta_2 = 0$ . The models $\beta_0 + \beta_1 x_1$ and $\beta_0 + \beta_2 x_2$ with two different predictors are in general not nested in each other, and a model for $Y$ cannot be compared in this way with a model for $\log(Y/y_0)$ , whose residuals are on a different scale. If adding a predictor with missing values silently drops records, the two fits no longer use the same observations and their values of $\text{SSE}$ cannot be compared.

Partial F test

Definition Let model 1 have $p_1$ coefficients, counting the intercept, and let model 0 be model 1 with $q$ of those coefficients set to zero, so it has $p_0 = p_1 - q$ . Fit both by least squares to the same $n$ observations of the same response, with design matrices of full column rank and $n > p_1$ . To test $H_0$ , that the $q$ omitted coefficients are all zero, against $H_1$ , that at least one of them is not, compute

F = \frac{(\text{SSE}_0 - \text{SSE}_1)/q}{\text{SSE}_1/(n - p_1)} .

Suppose the mean function of model 1 is correct and, given the predictors, the errors are independent and normal with mean $0$ and a common variance. Then under $H_0$ the statistic has the $F_{q,\, n - p_1}$ distribution, large values count against $H_0$ , and the p-value is the probability that an $F_{q,\, n - p_1}$ variable is at least the observed $F$ . When $q = 1$ the statistic equals the square of the $t$ statistic for the added coefficient in model 1, and the F test agrees with the two-sided $t$ test.

In plain words The numerator is the drop in squared error per added coefficient, and the denominator estimates the error variance from the larger model. Under $H_0$ each added term is expected to remove about one error variance of squared error, so values near $1$ are ordinary. The word partial means the extra terms are judged with the smaller model’s terms kept in place. Testing a block of terms jointly answers a different question from scanning their separate p-values and keeping the smallest one.

Example Twelve observations have $\text{SST} = 100$ . Model A has $2$ coefficients and $\text{SSE} = 40$ , and model B adds one coefficient and has $\text{SSE} = 30$ . Then

F = \frac{(40 - 30)/1}{30/9} = 3 .

If the classical conditions hold, the p-value with $1$ and $9$ degrees of freedom is $P(F_{1,9} \ge 3) \approx 0.117$ . At a significance level of $0.05$ chosen in advance, $H_0$ is not rejected. That does not show that the added coefficient is zero or that model A is true. It says only that these twelve observations do not give clear evidence of the added term’s contribution.

Adjusted R-squared

Definition Consider a least-squares model with an intercept and $p$ coefficients in all, fitted to $n > p$ observations with $\text{SST} > 0$ . Let $s^2 = \text{SSE}/(n - p)$ be its estimate of the error variance and $s_y^2 = \text{SST}/(n - 1)$ the sample variance of the response. The adjusted R-squared is

\bar{R}^2 = 1 - \frac{s^2}{s_y^2} ,

or equivalently

\bar{R}^2 = 1 - (1 - R^2)\,\frac{n - 1}{n - p} .

It never exceeds $R^2$ , and it is negative when $s^2 > s_y^2$ . For two nested models fitted to the same observations, $\bar{R}^2$ is higher for the larger model exactly when the partial F statistic comparing them is greater than $1$ , and the two values are equal when $F = 1$ .

In plain words Ordinary $R^2$ compares sums of squares, which never favor the smaller of two nested models. Adjusted $R^2$ compares variance estimates, dividing each sum of squares by its degrees of freedom, so an added coefficient has to reduce $\text{SSE}$ by enough to make up for the residual degree of freedom it uses. It describes fit to the data at hand with a penalty for size. It does not estimate the error on new data and is not a general rule for choosing a model.

Example For the twelve observations, $s_y^2 = 100/11$ . Model C adds one more coefficient to model B and has $\text{SSE} = 29$ .

Model	$p$	$\text{SSE}$	$R^2$	$\bar{R}^2$
A	$2$	$40$	$0.60$	$0.56$
B	$3$	$30$	$0.70$	$\approx 0.6333$
C	$4$	$29$	$0.71$	$0.60125$

Going from B to C raises $R^2$ but lowers $\bar{R}^2$ , in line with the partial F statistic $(30 - 29)/(29/8) \approx 0.276$ , which is below $1$ .

Multicollinearity

Definition Multicollinearity means that a predictor is close to a linear combination of the other predictors in the data. For a numerical predictor $x_j$ in a least-squares model with an intercept, regress $x_j$ on the other predictors, with an intercept, and call the $R^2$ of this auxiliary regression $R_j^2$ . When $R_j^2 < 1$ the variance inflation factor is

\operatorname{VIF}_j = \frac{1}{1 - R_j^2} .

If the errors are uncorrelated with common variance $\sigma^2$ , the estimate $b_j$ of the coefficient of $x_j$ has, given the predictor values,

\operatorname{Var}(b_j) = \frac{\sigma^2}{S_{jj}} \operatorname{VIF}_j ,

where $S_{jj}$ is the sum of squared deviations of $x_j$ from its mean. So $\operatorname{VIF}_j$ is the factor by which this variance exceeds its value when $x_j$ is uncorrelated with the other predictors, with $\sigma^2$ and $S_{jj}$ unchanged, and $\sqrt{\operatorname{VIF}_j}$ is the factor for the standard error.

In plain words When predictors move together in the data, the data say a lot about their combined association with the response and little about how to divide it among them. Individual coefficients then have inflated standard errors, and their values and even their signs can change a lot with small changes in the data, while fitted values for profiles like those in the data can stay stable. A VIF measures uncertainty, not bias. A coefficient with VIF $10$ is not ten times too large, and a rule-of-thumb cutoff is no substitute for asking why the predictors are related.

Example If $R_j^2 = 0.9$ , then $\operatorname{VIF}_j = 1/0.1 = 10$ , and the standard error of $b_j$ is $\sqrt{10} \approx 3.16$ times what it would be with $R_j^2 = 0$ , other things equal. With $R_j^2 = 0.5$ the factor on the variance is $2$ and on the standard error $\sqrt{2} \approx 1.41$ .

Identifiability

Definition A parameter of a model is identifiable when different values of it always give different distributions of the observable data. In a linear regression with design matrix $X$ the data depend on the coefficient vector $\beta$ only through the mean vector $X\beta$ and the error distribution, so the coefficients are identifiable exactly when the columns of $X$ are linearly independent. Under perfect collinearity some column is an exact linear combination of the others, distinct coefficient vectors give the same fitted values, and least squares has infinitely many solutions. Some combinations of the coefficients can still be identifiable.

In plain words If two settings of the coefficients make exactly the same predictions for every possible record of this kind, no amount of such data can tell them apart. Software usually drops a column or reports a missing estimate. More data with the same exact relation does not help. Data in which the relation fails, or a restriction chosen in advance, can.

Example Suppose $x_2 = 2x_1$ in every record of a model with an intercept. Then

\begin{aligned} &\beta_1 x_1 + \beta_2 x_2 \\[0.4em] &\quad = (\beta_1 + 2\beta_2)\, x_1 , \end{aligned}

so only the combination $\beta_1 + 2\beta_2$ enters the fitted values. The pairs $(\beta_1, \beta_2) = (2, 1)$ and $(4, 0)$ both give $4x_1$ , which at $x_1 = 0, 1, 2, 3$ equals $0$ , $4$ , $8$ and $12$ either way. The data identify $\beta_1 + 2\beta_2$ but not $\beta_1$ and $\beta_2$ separately.

Model selection

Definition Model selection is the use of data to choose among candidate models, for example by deciding which predictors, interactions or curved terms to include. The selection procedure is everything that was tried, such as tests, criteria and plots, together with the rule for picking a model. Standard p-values and confidence intervals assume a model fixed before the data were seen. When the model is chosen by searching the same data, they generally lose their stated properties, because they ignore the uncertainty of the selection itself.

In plain words Start from what the model is for. For prediction, judge the whole procedure, search included, on observations not used to fit or choose the model. For explaining an association, choose adjustment variables from knowledge of the design and the subject. Dropping a confounder because its p-value is large changes what the remaining coefficients mean. Record what was searched and label findings of the search as exploratory unless new data or a method that accounts for the selection confirms them, as with multiple testing.

Example Suppose $20$ candidate predictors are unrelated to the response and each is tested at the $0.05$ level. If the $20$ tests were independent, the probability that at least one comes out significant would be

1 - 0.95^{20} \approx 0.64 .

Reporting only the winner, with its p-value below $0.05$ , makes the result of a search look like a confirmed finding.

Check yourself

In a regression with an intercept, regressing the predictor $x_3$ on the other predictors gives $R^2 = 0.75$ . Find $\operatorname{VIF}_3$ and the factor by which the standard error of $b_3$ exceeds its value for an uncorrelated predictor. A colleague concludes that $b_3$ is four times too large. Is that right?

Show the answer

The variance inflation factor is $1/(1 - 0.75) = 4$ , so the sampling variance of $b_3$ is four times, and its standard error $\sqrt{4} = 2$ times, what it would be if $x_3$ were uncorrelated with the other predictors, with the same error variance and the same spread of $x_3$ . The colleague is wrong. The VIF describes how much $b_3$ would vary from sample to sample, not whether it is too large. If the linear model is correct, $b_3$ is still unbiased, but for a given estimate of the error variance a confidence interval for $\beta_3$ is twice as wide as it would be in the uncorrelated comparison.

For $n = 11$ observations with $\text{SST} = 100$ , a model with $3$ coefficients has $\text{SSE} = 40$ , and a model that adds one coefficient to it has $\text{SSE} = 35$ . Find $R^2$ and $\bar{R}^2$ for both models and the partial F statistic. Why do the adjusted values agree?

Show the answer

The values of $R^2$ are $1 - 40/100 = 0.60$ and $1 - 35/100 = 0.65$ . With $s_y^2 = 100/10 = 10$ , the adjusted values are $1 - (40/8)/10 = 0.5$ and $1 - (35/7)/10 = 0.5$ . The partial F statistic is

F = \frac{(40 - 35)/1}{35/7} = 1 .

The drop in $\text{SSE}$ , which is $5$ , exactly equals the larger model’s residual mean square, also $5$ , and $F = 1$ is precisely the case in which adjusted $R^2$ neither rises nor falls. Under the classical assumptions, $P(F_{1,7} \ge 1) \approx 0.351$ , so these data give no clear evidence for the added term, which is not the same as evidence that its coefficient is zero.

A model with an intercept uses predictors $x_1$ and $x_2$ , and every record has $x_2 = x_1 + 5$ . Since $x_2$ is not a multiple of $x_1$ , someone expects all three coefficients to be identifiable. Are they? Which combinations are identified?

Show the answer

They are not. Substituting gives

\begin{aligned} &\beta_0 + \beta_1 x_1 + \beta_2 (x_1 + 5) \\[0.4em] &\quad = (\beta_0 + 5\beta_2) \\[0.4em] &\qquad + (\beta_1 + \beta_2)\, x_1 . \end{aligned}

The columns of the design matrix are linearly dependent, since $x_2 - x_1 - 5 \cdot 1 = 0$ in every record. For example, the coefficient vectors $(\beta_0, \beta_1, \beta_2) = (1, 2, 1)$ and $(6, 3, 0)$ both give $6 + 3x_1$ . When $x_1$ takes at least two different values, the identified quantities are $\beta_0 + 5\beta_2$ , $\beta_1 + \beta_2$ and their linear combinations. The condition concerns linear dependence among all the columns, the intercept column included, not only whether one predictor is a multiple of another.

Comparing Regression Models and Multicollinearity

Key terms

Nested models

Partial F test

Adjusted R-squared

Multicollinearity

Identifiability

Model selection

Check yourself

How this lesson connects

Builds on