Lecture 12: Comparison of Multiple Linear Regression Models

Download the R markdown file for this lecture.

We have now considered two types of testing problems in multiple linear regression:

Testing whether the response is related to at least one explanatory variable;
Testing the effect of a given (single) explanatory variable having adjusted for other variables.

In this lecture, we investigate the importance of a group of covariates simultaneously.

This problem is equivalent to comparing a given model with a simplified version of that model.

Nested Models

A linear model M0 is said to be nested within another model M1 if M0 can be recovered as a special case of M1 by setting the parameters of M1 to the necessary values.

For the paramo data, define model M1 by

\[E[\mbox{N}] = \beta_0 + \beta_1 \mbox{AR} + \beta_2 \mbox{EL} + \beta_3 \mbox{DEc} + \beta_4 \mbox{DNI}\]

Then the model M0 defined by

\[E[\mbox{N}] = \beta_0 + \beta_1 \mbox{AR} + \beta_3 \mbox{DEc}\] is nested within M1 because we can obtain M0 by setting $\beta_2 = \beta_4 = 0$ in M1.

F Tests for Nested Models

The general ideas about model comparison from lecture 11 continue to apply.
- We want a “cheap” model (i.e. one with few parameters);
- We want a model that fits well (i.e. one with small RSS.
The choice of model will be a question of whether improvement in goodness of fit of complex model M1 over simpler model M0 is worth the cost.

Selecting a model is equivalent to testing hypotheses about parameters of M1 (more complex model).

Suppose that linear model M0 has q explanatory variables (hence q+1 regression parameters, including intercept).
Suppose that M0 is nested within linear model M1 which has p > q explanatory variables (hence p+1 regression parameters, including intercept).

Comparison of the models can be achieved by testing

H₀: $\beta_j = 0$ for all $j \in J$ (i.e. M0 adequate); versus

H₁: $\beta_j$ not zero for all $j \in J$ (i.e. M1 better)

where J indexes the p-q variables that appear in M1 but not M0.

The F-statistic to test these hypotheses is

\[F = \frac{[RSS_{M0} - RSS_{M1}]/(p-q)}{RSS_{M1}/(n-p-1)}\]

As before, extreme, large values of F provide evidence against H₀ (hence evidence that we should prefer model M1 to M0).

If H₀ is correct then F has an F distribution with (p-q),(n-p-1) degrees of freedom.

Hence if f_obs is the observed value of the test statistic, and X is a random variable from an F_p-q,n-p-1 distribution, then the P-value is given by \[P= P(X \ge f_{obs})\]

Interpretation of Test Results

As we have seen earlier, the importance of variables in a multiple regression is dependent upon context (i.e. what other variables are taken into account).
Retention of H₀ (i.e. acceptance that model M0 is adequate) does not mean that the p-q variables indexed by J are unrelated to the response.
Rather, retention of H₀ indicates that the variables indexed by J do not provide additional information about the response having adjusted for all the other (q) variables.

Model Comparison for the Paramo Data

Suppose (for the purposes of argument) that a scientist theorizes that EL and DNI provide no further information about N once we have adjusted for AR and DEc.
Define M1 by \[E[\mbox{N}] = \beta_0 + \beta_1 \mbox{AR} + \beta_2 \mbox{EL} + \beta_3 \mbox{DEc} + \beta_4 \mbox{DNI}\] and M0 by \[E[\mbox{N}] = \beta_0 + \beta_1 \mbox{AR} + \beta_3 \mbox{DEc}\] as in the previous example.
Testing the scientist’s theory requires that we test H₀: \[\beta_2 = \beta_4 = 0~~~~\mbox{versus}~~~~~H_1: \beta_2, \beta_4 \mbox{ not both zero.}\]

Download paramo.csv

## Paramo <- read.csv(file = "https://r-resources.massey.ac.nz/161221/data/paramo.csv", 
##     header = TRUE, row.names = 1)
Paramo.lm0 <- lm(N ~ AR + DEc, data = Paramo)
anova(Paramo.lm0)

Analysis of Variance Table

Response: N
          Df Sum Sq Mean Sq F value   Pr(>F)   
AR         1 508.92  508.92  12.937 0.004193 **
DEc        1 557.23  557.23  14.165 0.003134 **
Residuals 11 432.71   39.34                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Paramo.lm1 <- lm(N ~ AR + EL + DEc + DNI, data = Paramo)
anova(Paramo.lm1)

Analysis of Variance Table

Response: N
          Df Sum Sq Mean Sq F value   Pr(>F)   
AR         1 508.92  508.92 11.3208 0.008328 **
EL         1  45.90   45.90  1.0211 0.338661   
DEc        1 537.39  537.39 11.9541 0.007189 **
DNI        1   2.06    2.06  0.0457 0.835412   
Residuals  9 404.59   44.95                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

For model M0, calculations give RSS_M0 = 432.7 and q=2.
For model M1, calculations give RSS_M1 = 404.6 and p=4.
Hence F-statistic is \[F = \frac{[RSS_{M0} - RSS_{M1}]/(p-q)}{RSS_{M1}/(n-p-1)} = \frac{[432.7 - 404.6]/(4-2)}{404.6/(9)} = 0.31\]
Corresponding P-value is right hand tail probability: \[P(X \ge 0.31) = 0.74~~~~~~~~~~\mbox{where $X \sim F_{2,9}$}\]
We conclude that H₀ can be retained (i.e. model M0 is adequate). There is no evidence that EL and DNI provide further information about N once we have adjusted for AR and DEc.

Model Comparison Using R

anova(Paramo.lm0, Paramo.lm1)

Analysis of Variance Table

Model 1: N ~ AR + DEc
Model 2: N ~ AR + EL + DEc + DNI
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     11 432.71                           
2      9 404.59  2     28.12 0.3128 0.7391

Connection with the Omnibus F Test

The omnibus F test of model fit discussed in lecture 10 is a particular case of the more general methodology for model comparison presented in this lecture.
Specifically, in omnibus F test we always take M0 to be the null model with q=0 explanatory variables, and M1 to be the full model with p explanatory variables.

Connection with T Tests for Single Variables

F tests can be used to compare two models that differ by a single explanatory variable.
T test can also be used to assess the importance of any given explanatory variable.
There is a connection between the two approaches. Suppose we test H₀: \[\beta_j = 0~~~~\mbox{versus}~~~~H_1:\beta_j \ne 0\] having adjusted for certain other variables. If the observed t-test statistic is t_obs, and the observed F test statistic is f_obs, then f_obs = t_obs².
The P-value from the two tests will be the same.

More Testing for the Paramo Data

Suppose that we wish to test for the importance of EL having adjusted for AR (only).
Consider model M1 given by \[E[\mbox{N}] = \beta_0 + \beta_1 \mbox{AR} + \beta_2 \mbox{EL}\] and M0 by \[E[\mbox{N}] = \beta_0 + \beta_1 \mbox{AR}\]
We want to test H₀: \[\beta_2 = 0~~~~\mbox{versus}~~~~H_1:\beta_2 \ne 0\] which is equivalent to comparing models M0 and M1.
Will perform t and F tests using R.

Paramo.lm2 <- lm(N ~ AR, data = Paramo)
Paramo.lm3 <- lm(N ~ AR + EL, data = Paramo)
anova(Paramo.lm2, Paramo.lm3)

Analysis of Variance Table

Model 1: N ~ AR
Model 2: N ~ AR + EL
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     12 989.94                           
2     11 944.04  1    45.901 0.5348 0.4799

The F-statistic (from the anova() function output) is f = 0.5348 with corresponding P-value P=0.4799.
The t-statistic (from the summary() function output) is t=0.731 with corresponding P-value P=0.4799.
Note that f = 0.5348 = 0.731² = t².